{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Difference-in-Difference\n",
    "\n",
    "## Three Billboards in the South of Brazil\n",
    "\n",
    "I remember when I worked with marketing and a great way to do it was with internet advertisement. Not because it is very efficient (although it is), but because it is very easy to know if it's effective or not. With online marketing, you have a way of knowing which customers saw the ad and you can track them with cookies to see if they ended up on your landing page or clicked some download button. You can also use machine learning to find prospects that are very similar to your customers and present the ad only to them. In this sense, online marketing is very precise: you target only those you want to and you can see if they respond as you would like them to. \n",
    "\n",
    "But not everyone is susceptible to online marketing. Sometimes you have to resort to less precise techniques, like a TV campaign or placing a billboard down the street. Usually, diversity of marketing channels is something marketing departments look for. But if online marketing is a professional fishing rod to catch that specific type of tuna, billboard and TV are giant nets you throw at a fish shoal and hope to catch at least some big ones. Another problem with billboard and TV ads is that it is much harder to know how effective they are. Sure, you could measure the purchase volume, or whatever you want to drive, before and after placing a billboard somewhere. If there is an increase, there is some evidence that the marketing is effective. But how would you know if this increase is not just some natural trend in the awareness of your product? In other words, how would you know the counterfactual \\\\(Y_0\\\\) of what would have happened if you didn't set up the billboards in the first place? \n",
    "\n",
    "![img](./data/img/diff-in-diff/secrets.png)\n",
    "\n",
    "One technique to answer these types of questions is simple Difference-in-Difference, or diff-in-diff for close friends. Diff-in-diff is commonly used to assess the effect of macro interventions, like the effect of immigrants on unemployment, the effect of gun law changes in crime rates or simply the difference in user engagement due to a marketing campaign. In all these cases, you have a period before and after the intervention and you wish to untangle the impact of the intervention from a general trend. As a motivating example, let's look at a question similar to the one I had to answer.\n",
    "\n",
    "In order to figure out how good billboards were as a marketing channel, we placed 3 billboards in the city of Porto Alegre, the capital of the state of Rio Grande do Sul. We wanted to see if that boosted deposits into our savings account. As a note for those not very familiar with Brazilian geography, Rio Grande do Sul is part of the south of the country, one one of the most developed regions. \n",
    "\n",
    "Having this in mind, we decided to also look at data from another capital from the south, Florianopolis, the capital city of the state of Santa Catarina. The idea is that we could use Florianopolis as a control sample to estimate the counterfactual \\\\(Y_0\\\\) when compared to Porto Alegre. (By the way, this was not the true experiment, which is confidential, but the idea is very similar). We placed the billboard in Porto Alegre for the entire month of June. The data we have looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import style\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "style.use(\"fivethirtyeight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>deposits</th>\n",
       "      <th>poa</th>\n",
       "      <th>jul</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>42</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>52</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>119</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>21</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   deposits  poa  jul\n",
       "0        42    1    0\n",
       "1         0    1    0\n",
       "2        52    1    0\n",
       "3       119    1    0\n",
       "4        21    1    0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv(\"data/billboard_impact.csv\")\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that deposits are our outcome variable, the one we wish to increase with the billboards. POA is a dummy indicator for the city of Porto Alegre. When it is zero, it means the samples are from Florianopolis. Jul is a dummy for the month of July, or for the post intervention period. When it is zero it refers to samples from May, the pre-intervention period.\n",
    "\n",
    "## DID Estimator\n",
    "\n",
    "To avoid confusion between Time and Treatment, from now on, I'll use D to denote treatment and T to denote time. Let \\\\(Y_D(T)\\\\) be the potential outcome for treatment D on period T. In an ideal world where we have the ability to observe the counterfactual, we would estimate the treatment effect of an intervention the following way:\n",
    "\n",
    "$\n",
    "\\hat{ATET} = E[Y_1(1) - Y_0(1)|D=1]\n",
    "$\n",
    "\n",
    "In words, the causal effect is the outcome in the period post intervention in case of a treatment minus the outcome in also in the period after the intervention, but in the case of no treatment. Of course, we can't measure this because \\\\(Y_0(1)\\\\) is counterfactual. \n",
    "\n",
    "One way around this is a before and after comparison.\n",
    "\n",
    "$\n",
    "\\hat{ATET} = E[Y(1)|D=1] - E[Y(0)|D=1]\n",
    "$\n",
    "\n",
    "In our example, we would compare the average deposits from POA before and after the billboard was placed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "41.04775"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "poa_before = data.query(\"poa==1 & jul==0\")[\"deposits\"].mean()\n",
    "\n",
    "poa_after = data.query(\"poa==1 & jul==1\")[\"deposits\"].mean()\n",
    "\n",
    "poa_after - poa_before"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This estimator is telling us that we should expect deposits to increase R$ 41,04 after the intervention. But can we trust this?\n",
    "\n",
    "Notice that \\\\(E[Y(0)|D=1]=E[Y_0(0)|D=1]\\\\), that is, the observed outcome for the treated unit **before the intervention** is equal to the counterfactual outcome for the treated unit also before the intervention. Since we are using, this to estimate the counterfactual **after the intervention**  \\\\(E[Y_0(1)|D=1]\\\\), this estimation above assumes that \\\\(E[Y_0(1)|D=1] = E[Y_0(0)|D=1]\\\\). \n",
    "\n",
    "It is saying that in the case of no intervention, the outcome in the latter period would be the same as the outcome from the starting period. This would obviously be false if your outcome variable follows any kind of trend. For example, if deposits are going up in POA, \\\\(E[Y_0(1)|D=1] > E[Y_0(0)|D=1]\\\\), i.e. the outcome of the latter period would be greater than that of the starting period even in the absence of the intervention. With a similar argument, if the trend in Y is going down, \\\\(E[Y_0(1)|D=1] < E[Y_0(0)|D=1]\\\\). This is to show that this before and after thing is not a great estimator. \n",
    "\n",
    "Another idea is to compare the treated group with an untreated group that didn't get the intervention:\n",
    "\n",
    "$\n",
    "\\hat{ATET} = E[Y(1)|D=1] - E[Y(1)|D=0]\n",
    "$\n",
    "\n",
    "In our example, it would be to compare the deposits from POA to that of Florianopolis in the post intervention period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-119.10175000000001"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fl_after = data.query(\"poa==0 & jul==1\")[\"deposits\"].mean()\n",
    "poa_after - fl_after"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This estimator is telling us that the campaign is detrimental and that customers will decrease deposits by R$ 119.10. \n",
    "\n",
    "Notice that \\\\(E[Y(1)|D=0]=E[Y_0(1)|D=0]\\\\). And since we are using \\\\(E[Y(1)|D=0]\\\\) to estimate the counterfactual for the treated after the intervention, we are assuming we can replace the missing counterfactual like this: \\\\(E[Y_0(1)|D=0] = E[Y_0(1)|D=1]\\\\). But notice that this would only be true if both groups have a very similar baseline level. For instance, if Florianopolis has way more deposits than Porto Alegre, this would not be true because \\\\(E[Y_0(1)|D=0] > E[Y_0(1)|D=1]\\\\). On the other hand, if the level of deposits are lower in Florianopolis, we would have \\\\(E[Y_0(1)|D=0] < E[Y_0(1)|D=1]\\\\). \n",
    "\n",
    "Again, this is not a great idea. To solve this, we can use both space and time comparison. This is the idea of the difference in difference approach. It works by replacing the missing counterfactual the following way:\n",
    "\n",
    "$\n",
    "E[Y_0(1)|D=1] = E[Y_1(0)|D=1] + (E[Y_0(1)|D=0] - E[Y_0(0)|D=0])\n",
    "$\n",
    "\n",
    "What this does is take the treated unit before the treated unit **before the intervention** and adds a trend component to it, which is estimated using the control \\\\(E[Y_0(1)|T=0] - E[Y_0(0)|T=0]\\\\). In words, it is saying that the treated **after the intervention**, had it not been treated, would look like the **treated before the treatment** plus a growth factor that is the same as the growth of the control. \n",
    "\n",
    "It is important to notice that this assumes that the trends in the treatment and control are the same:\n",
    "\n",
    "$\n",
    "E[Y_0(1) − Y_0(0)|D=1] = E[Y_0(1) − Y_0(0)|D=0]\n",
    "$\n",
    "\n",
    "where the left hand side is the counterfactual trend. Now, we can replace the estimated counterfactual in the treatment effect definition \\\\(E[Y_1(1)|D=1] - E[Y_0(1)|D=1]\\\\)\n",
    "\n",
    "$\n",
    "\\hat{ATET} = E[Y(1)|D=1] - (E[Y(0)|D=1] + (E[Y(1)|D=0] - E[Y(0)|D=0])\n",
    "$\n",
    "\n",
    "If we rearrange the terms, we get the classical Diff-in-Diff estimator.\n",
    "\n",
    "$\n",
    "\\hat{ATET} = (E[Y(1)|D=1] - E[Y(1)|D=0]) - (E[Y(0)|D=1] - E[Y(0)|D=0])\n",
    "$\n",
    "\n",
    "It gets that name because it gets the difference between the difference between treatment and control after and before the treatment. \n",
    "\n",
    "Here is what that looks in code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.524557692307688"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fl_before = data.query(\"poa==0 & jul==0\")[\"deposits\"].mean()\n",
    "\n",
    "diff_in_diff = (poa_after-poa_before)-(fl_after-fl_before)\n",
    "diff_in_diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Diff-in-Diff is telling us that we should expect deposits to increase by R$ 6.52 per customer. Notice that the assumption that diff-in-diff makes is much more plausible than the other 2 estimators. It just assumes that the growth pattern between the 2 cities are the same. But it doesn't require them to have the same base level nor does it require the trend to be zero. \n",
    "\n",
    "To visualize what diff-in-diff is doing, we can project the growth trend from the untreated into the treated to see the counterfactual, that is, the number of deposits we should expect if there were no intervention."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApcAAAErCAYAAACLjKgCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3hkVZ3v//d31zXp3NNJ03Y3NEq4i42jyE1AQAVvjOPlACqIeBnhKKPDDDI6oyMiyPiTUTnjFQSkf4gKas8chnNwFIFxgFF0AGmYINduoJPu3Dup1GWv80ddUlWpJJVkJ6kkn9fz1LOrdu3atap5gE+vtb5rmXMOEREREZEgeEvdABERERFZORQuRURERCQwCpciIiIiEhiFSxEREREJjMKliIiIiAQmvBRfOjg4qBJ1ERERkWWuubnZys+p51JEREREAqNwKSIiIiKBWTXhsru7e6mbICIiIrIgainnrJpwKSIiIiILT+FSRERERAKzJNXiU3HOMTIygu/7gd87Ho8zODgY+H1XO8/zaGhowGxSsZiIiIisQjUVLkdGRojFYkSj0cDvHYvFiMfjgd93tUsmk4yMjNDY2LjUTREREZEaUFPD4r7vL0iwlIUTjUYXpKdZRERElqeaCpciIiIisrzV1LB4LWhra+PQQw8tvN66dSvPPvss11xzDbfccssStkxERESk9ilclqmrq+Pee+8tOffss88uUWtEREREpjaS8knX2Ow0hUsRERGRGpLMOHoTPj1jGXaNZegZ89k1mjuOlR73ph3nHVTPBR1L3eoJNR0uW763M9D7DZy3YcZrxsbGOP744wHYb7/92Lp1a6BtEBERkdXHd46+cZ9do/nQOPnYmzv2jVffFRkPge8WsOFzUNPhcilUGhYXERERKeecYzjl6J3Uozg5PPaM+WSqDIGeQWfco7MuxLo6j8763LFu4tiZOzZFDDOju3v3wv7YWajpcFlNT2O1EolEYPcSERGRlWs84wqBsHwYunx4eqzaxAi0xozOeDYYrqvPHetCJaFxXZ1HW8wj5C3fzUlqOlyKiIiIBCHjO/aM+xM9iqMZehOVQ+NAsvrAWBcy1tXnQ+LUx466ELHQ8g2Ms6FwWaVf/epXJUsUXX/99Rx11FFL2CIREZHVzTnHUMpNDD+PVp7L2DOWDZLVzk0MG4Vh5/Jh6JLQWO/REDZtgVxG4bLMzp2Ti4he+9rX8uKLLy5Ba0RERFafsfRUw9LZ58WhMZGp/r5tMW/S3MVKcxpbYx6eAuOczRguzWwTcCOwD+AD33bOfdXM2oBbgM3A08C7nXP9lo3vXwXeBIwC73fOPbgwzRcREZHlIOM7dicqzGGscByaxbB0Q9hK5jBO6l3MneuIe0RXybD0Uqum5zIN/KVz7kEzawR+a2Z3Au8H/s05d6WZfQr4FHAJcDrQlXu8BvhG7igiIiIriHOOwaRjV/kw9GimJCz2Jnx2z2JYOuKRLXypr9DLmDuuqwvRUefRENFO1rVmxnDpnHsBeCH3fNjMtgMbgDOAk3KX3QDcRTZcngHc6JxzwH1m1mJm63P3ERERkRo3mvYnDT9PNacxWeWSjAasjXvTFL2ECoUxLVHNY1zOZjXn0sw2A0cC9wPr8oHROfeCmXXmLtsAPFf0sR25cxXDZXd3d+F5PB4nFovNpkmzouWIFsbQ0BA9PT1L3QwREZlG2kF/0tiTgj1Jyz5SZcfc872Z6oPdmpCjPepoj0xxjDraI9AacYSn6mT0gb3Zx26yD5m94ky1kLq6uqZ9v+pwaWYNwK3AXzjnhqb5G0WlN6bsCC9u4ODgIPF4vNomzUoikViwe692TU1NbNq0aambISKy6jjn6C9eXic/DF1hLuOehD/1/4zLRD1mLHrJV1DXT5kYZTF1d3fPGPoWS1Xh0swiZIPlVufcbbnTu/LD3Wa2Hsh3Xe0AipPGRuD5oBosIiKy0u1N+YVgONXSOj2jPj2JDKlZDEuXL69TfOwomsvYrGFpmYdqqsUNuBbY7pz7StFb24BzgStzx58Vnf+fZvYDsoU8g8tpvmVbWxuHHnoomUyGAw88kG984xvU19ezc+dOLr74Yh5//HF83+eNb3wjl112GdFotPDZs846i927d3PnnXcu4S8QEZFalMw4ehP5gFh5t5f8Ujsj6eqrpZuiNuMC3uvqQrTHPcLLeNcXWT6q6bk8Dngf8LCZ/T537m/Ihsofmtn5wLPAu3Lv3U52GaInyC5FdF6gLV5gxXuLf+hDH+K6667jwgsv5H3vex8f+MAHuPnmm8lkMlx00UVcdtllXHbZZQAMDAzw0EMPsWbNGp5++mk2b968hL9CREQWg19hWHqqhbz7xqvsYgTiofJh6bKwmF92Jx4iHlZglNpSTbX4vVSeRwlwSoXrHXDhPNtVE4455hj+8Ic/cPfddxOLxXjve98LQCgU4otf/CKveMUruPTSS6mvr2fbtm2cdtppdHR0cNttt/HJT35yiVsvIiJz4ZxjJO3oGZ28gPekXV/GfKrtZPQMOuOTh6ErzWlsimhYWpavmt6hp+Hck4K7FzByw11VX59Op7nzzjs59dRT2b59O1u2bCl5v6mpiY0bN/Lkk09y+OGHc+utt3LJJZfQ2dnJOeeco3ApIlJjxjOlu75UKnrJH0dnMSzdUjwsXV+2tE5ReGyLeYQ0LC2rQE2Hy6UwNjbG8ccfD2R7Lt/3vvdx7bXXVvwbpHMOM6Onp4cnn3ySY445BjMjHA7z6KOPluxFLiIiwcv4jr7yYekp5jQOzGLXl7qQFdZcnGoOY0cuOMa064tIiZoOl7PpaZxJIpGgmoWIiudc5h1yyCFs27at5NzQ0BA7d+5k//3358Ybb2RwcJAjjjgCgOHhYW677TaFSxGROXDOMZRyU85h7C0elk74ZKrMjGHLVktPGo4un9NY79EQ1rC0yFzVdLisFSeeeCKf+9znuPnmmznrrLPIZDJ85jOf4eyzz6a+vp5bb72VH//4xxx11FEAPP3007z97W/nM5/5zBK3XESkdiTSjp5E5Srp8jmNiUz1922LeRWX1ilfl7E15uEpMIosOIXLKpgZN910E3/5l3/JP/zDP+D7Pq9//ev5u7/7O5555hl27NjBq1/96sL1mzdvprGxkd/85je86lWvWsKWi4gsrIzv2J2Yeu5i8XFoFsPSDWErrYquC9EZnzynsSPuEdWwtEhNsWxx9+IaHBys+KWDg4M0NzcvyHdqh56Fs5D/3ERk8TnnGEy6itXR5cfdCR9/FsPS+WHnqRbyzs9lbIho1xeR2ViqHXqam5sn/e1OPZciIqvEaNovLNQ907qMyeqXZGRt3JtmAe/8c48WDUuLrAoKlyIiy1jaz+76UjyHsfx1PlAOpWax60vECtXQU1ZM14dYG/eIaHkdESmicCkiUmOccwwkXXZf6dGypXXKjnsSPtVGxqhXvuvL5KKXfE9jfVjD0iIyNwqXIiKLZG/KLwTDSnMYe4uCY6rKYWmDQsHLuqLCl/LQuK4uRHNUy+uIyMJTuBQRmYeU7+gtCojlS+sUL7szMotdX5qKd32Z5tge9whrWFpEaojCpYhIGd85+st3falQ9LJrzKdvvPrKl1goWy1dceHuoiHpzroQdWEFRhFZnhQuy+zatYtLL72UBx98kFgsxr777ssVV1zBAQccEMj977nnHqLRKK95zWtm/dnzzz+f7du38573vIcLL7yw6s8NDAzw4x//mA9+8IOz/s68j370o5x22mmcccYZc76HyFJyzjGSdvSMlha5VFrAu2fMp9pORs+gIz710jrFw9NNEQ1Li8jKp3BZxDnHe9/7Xs466yyuu+46AB566CF6enoCC5f33nsvDQ0NswqX6XSaPXv2cP/99/PII4/M+jsHBwe59tpr5xUuRWpVMuPKQuLUx9FZDEu3FA9Lly3cXRwe22IeIQ1Li4gUKFwWufvuuwmHw3zgAx8onDviiCNwzvG3f/u3/PznP8fMuPjii/mzP/sz7rnnHq655hpuueUWAP7qr/6KLVu28J73vIeXv/zlnHXWWdxxxx2k02muv/56YrEY3/ve9wiFQtxyyy1cddVVHHjggXziE59gx44dAFxxxRUcffTRXHHFFbz44os8++yztLe38+ijj7J7926OP/54rrrqKrq7u7n++utJJpO89KUv5Vvf+hb19fX09PTwiU98gqeffhqAr3zlK3zrW9/iqaee4vjjj+d1r3sdb3jDG6Zs95e+9CXuuOMOEokERx11FP/4j/+onhZZdL5z7En4k4ahJ6/RmKF/vPrAWBcy1tVPHoauNDwd064vIiJzUtPhcu8vTpvV9V7jAdS9+ppJn19z8h1VfX779u1s2bJl0vlt27bx8MMPc++997Jnzx5OPvlkjj322Bnv197ezt133813v/tdvv71r/P1r3+d8847j4aGBj72sY8B8MEPfpALLriAY445hueee453vOMdPPDAAwD8/ve/54477qCuro5nnnmGM888k3vvvReAgw8+mHPPPReAL3zhC3z/+9/nIx/5CJdccgnHHXccW7duJZPJMDIywmc/+1m2b99e+Ow999wzZZs//OEPc8kllxSe33HHHZx++ulV/fmJTMc5x3Bqml1fiuY09iZ8MlVmxpBNrpauOKex3qMhrGFpEZGFVtPhslbcd999vOMd7yAUCtHZ2cmxxx7Lgw8+SGNj47Sfe+tb3wrAli1b+Od//ueK19x111089thjhdfDw8MMDw8DcPrpp1NXV1fxc48++iiXX345g4ODjIyMcMoppwDZ3tdvfvObAIRCIZqbmxkYGKj6t95999187WtfY2xsjP7+fg455BCFS5lWIu3oSZRWRVeqmO4Z8xmrNjECbTFvyjmMnfk5jvXZYWnt+iIiUjtqOlxW2+MY1OcPOeQQfvazn006P9X+6+FwGN+fqBRNJBIl78diMSAb8tLpdMV7+L7PnXfeWTFErlmzZsq2XnDBBWzdupWXv/zlbN26tdArWY2p2p1IJLj44ov55S9/ycaNG7niiism/SZZHTK+Y3fCpyeR61msGBqzx8Fk9YFxTdim3OmleE5jR9wjqmFpEZFlqabD5WI74YQT+PznP88NN9xQGHJ+8MEHaWlp4Sc/+Qlnn302/f39/PrXv+ayyy4jlUrx2GOPMT4+TiKR4Fe/+hVHH330tN/R0NBQ6JkEOPnkk/nOd77Dxz/+cSBbQHTEEUfM2NaRkRH22WcfUqkUP/rRj1i/fj0AJ554Itdeey0XXHABmUyGvXv30tjYWPKdmzZtqtjufJBsb29nZGSEbdu28ba3vW12f4hSs5xzDCanGZYuOu5O+PhVZsawURh2nqpiel1diI46j4aIdn0REVnpFC6LmBk33XQTl156KVdffTXxeLywFNHevXs5/vjjMTM+//nPs27dOgDe/va3c9xxx/Gyl72sqlB4+umnc84553D77bdz1VVX8aUvfYmLL76YY489lkwmw7HHHsvVV189430+/elPc8opp7Bp0yYOPfRQRkZGALjyyiu56KKLuOmmm/A8j6985SscddRRHH300RxzzDGceuqpXHbZZRXb3dLSwrnnnsuxxx7Lvvvuy5FHHjmPP01ZLGNpV3EpnUqhcTxT/X3Xxr0pi16Kg2OLhqVFRKSITTXku5AGBwcrfung4CDNzc0L8p2JRIJ4PL4g917tFvKf22qVzg1LTxqGLhue7hnLMJSq/t/hxogVFb9MsetLfYi1cY+IltcREVk2uru76erqWvTvbW5unvQ/C/VciiwS5xwDyWwv467RiaV0Ji/m7bM74VNtZIx6TLtwd2c8Gxg74h5rNCwtIiILbMZwaWbXAW8Bepxzh+fObQG+CcSBNHCBc+4By67x8VXgTcAo8H7n3IML1XiRWrA35dObyPYqTjWHMR8cU1XuFGjkd30p61msnzyXsTmq5XVERKR2VNNzeT1wDXBj0bmrgL93zv2rmb0p9/ok4HSgK/d4DfCN3FFkWUn5jt6igFhpaZ1do9n1GIdnMSzdVLzryzTH9rhHWMPSIiKyDM0YLp1zd5vZ5vLTQFPueTPwfO75GcCNLjuR8z4zazGz9c65FwJqr8ic+c7RP57d9aW3YmicWMh7z3iVXYxALERJkctUcxo760LUhRUYRURkZZvrnMu/AP6PmX0Z8ID8djUbgOeKrtuROzdluOzu7p5oTDiMmRGNRufYrOlpzcbgJZNJ+vr66OnpWbI2jGZgT9Kyj1TZMWnsSVE4l3HVhTsPR2sE2qMu+4iUHXPP10Yda0Iw5ah0MvtIDmb/ZRAREVkoxZlqIc1UODTXcPlR4BPOuVvN7N3AtcCpZKeKlZt2zLC4gc45RkZGGB8fn2OzpjY0NERTU9PMF8qseJ7H5s2bA5/zl8y4SUUuu8Yy9Fbobdybrn5YuqV4WLps4e7iLQTbYx4hDUuLiMgysVTV4pXMNVyeC1yUe/4j4Lu55zuATUXXbWRiyHxGZjbjlopz1dPTw6ZNm2a+UBaM7xx9435RpfRUxwz949UHxrqQTVpKp9Jcxs66EDHt+iIiIrKg5hounwdOBO4CTgby/bDbgP9pZj8gW8gzqPmWK5tzjuHUNLu+5OYw9ub2nq52a+mQUdKTWHosrp72aAirWlpERKRWVLMU0c1kK8HXmtkO4LPAh4CvmlkYSAAfzl1+O9lliJ4guxTReQvQZlkE41MMS1dayHus2sQItMW8CoUvZcd6jzbt+iIiIrIsVVMtftYUb/1JhWsdcOF8GyULI+M79uSqpXvGMpNCYvFxMFl9YKwPW2HNxamHp7OLeEc1LC0iIrKiaYeeZc45x2DSTVqwu9IwdW/Cx68yM4aNwrBzpd7F4hDZoF1fREREJEfhskaNpacalp4cGscz1d+3PT8sXaHopThAtmhYWkREROZA4XIRpX3H7sTUcxh7EhPnh2YxLN0YsSkX7i4eol4b94hoeR0RERFZQAqX85Qflt5VsUq6NETuTvjTL/pZJOpR6EXsqFgxnQ2MHXGPNRqWFhERkRqhcDmF0bRf6FWcaj3G/NzGZJU7BRrQEfcm9yzWl4XGuhDNUS2vIyIiIsvPqgqXKd/RWxQMK1VJ598bTlU/LN1UvOtL0bGj7PXauEdYw9IiIiKygq2acHn27+J031v1ZkHEQlReg7HCri91YQVGEREREVhF4TJsDs/yw9KV5zDmh6c74hqWFhEREZmLVRMu/+nwcY446ABCGpYWERERWTCrpsy4IYyCpYiIiMgCWzXhUkREREQWnsKliIiIiARG4VJEREREAqNwKSIiIiKBUbgUERERkcAoXIqIiIhIYBQuRURERCQwCpciIiIiEhiFSxEREREJjMKliIiIiARG4VJEREREAqNwKSIiIiKBmTFcmtl1ZtZjZo+Unf+YmT1uZn8ws6uKzl9qZk/k3nvjQjRaRERERGpTuIprrgeuAW7MnzCz1wFnAEc458bNrDN3/lDgTOAw4CXAz83sQOdcJuiGi4iIiEjtmbHn0jl3N9BXdvqjwJXOufHcNT2582cAP3DOjTvnngKeAI4KsL0iIiIiUsOq6bms5EDgtWZ2OZAALnbO/SewAbiv6LoduXNT6u7unmMTZm8xv0tERERkMS1Wzunq6pr2/bmGyzDQChwNvBr4oZm9FLAK17r5NDAo3d3di/ZdIiIiIouplnLOXKvFdwC3uawHAB9Ymzu/qei6jcDz82uiiIiIiCwXcw2XPwVOBjCzA4EosBvYBpxpZjEz2x/oAh4IoqEiIiIiUvtmHBY3s5uBk4C1ZrYD+CxwHXBdbnmiJHCuc84BfzCzHwKPAmngQlWKi4iIiKweM4ZL59xZU7z13imuvxy4fD6NEhEREZHlSTv0iIiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQmMwqWIiIiIBEbhUkREREQCo3ApIiIiIoFRuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQmMwqWIiIiIBEbhUkREREQCo3ApIiIiIoFRuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQmMwqWIiIiIBEbhUkREREQCM2O4NLPrzKzHzB6p8N7FZubMbG3utZnZ18zsCTN7yMxeuRCNFhEREZHaVE3P5fXAaeUnzWwT8Hrg2aLTpwNduceHgW/Mv4kiIiIislzMGC6dc3cDfRXeuhr4a8AVnTsDuNFl3Qe0mNn6QFoqIiIiIjUvPJcPmdnbgJ3Ouf8ys+K3NgDPFb3ekTv3wlT36u7unksT5mQxv0tERERkMS1Wzunq6pr2/VmHSzOrBz4NvKHS2xXOuQrnCmZqYFC6u7sX7btEREREFlMt5Zy59Fy+DNgfyPdabgQeNLOjyPZUbiq6diPw/HwbKSIiIiLLw6yXInLOPeyc63TObXbObSYbKF/pnHsR2Aack6saPxoYdM5NOSQuIiIiIitLNUsR3Qz8B3CQme0ws/Onufx24EngCeA7wAWBtFJEREREloUZh8Wdc2fN8P7moucOuHD+zRIRERGR5Ug79IiIiIhIYBQuRURERCQwCpciIiIiEhiFSxEREREJjMKliIiIiARG4VJEREREAqNwKSIiIiKBUbgUERERkcAoXIqIiIhIYBQuRURERCQwCpciIiIiEhiFSxEREREJjMKliIiIiARG4VJEREREAqNwKSIiIiKBUbgUERERkcAoXIqIiIhIYBQuRURERCQwCpciIiIiEhiFSxEREREJjMKliIiIiARmxnBpZteZWY+ZPVJ07h/M7DEze8jMfmJmLUXvXWpmT5jZ42b2xoVquIiIiIjUnmp6Lq8HTis7dydwuHPuCOC/gUsBzOxQ4EzgsNxn/snMQoG1VkRERERq2ozh0jl3N9BXdu7/OufSuZf3ARtzz88AfuCcG3fOPQU8ARwVYHtFREREpIYFMefyA8C/5p5vAJ4rem9H7pyIiIiIrALh+XzYzD4NpIGt+VMVLnPT3aO7u3s+TZiVxfwuERERkcW0WDmnq6tr2vfnHC7N7FzgLcApzrl8gNwBbCq6bCPw/HwaGJTu7u5F+y4RERGRxVRLOWdOw+JmdhpwCfA259xo0VvbgDPNLGZm+wNdwAPzb6aIiIiILAcz9lya2c3AScBaM9sBfJZsdXgMuNPMAO5zzv25c+4PZvZD4FGyw+UXOucyC9V4EREREaktM4ZL59xZFU5fO831lwOXz6dRIiIiIrI8aYceEREREQmMwqWIiIiIBEbhUkREREQCo3ApIiIiIoFRuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQmMwqWIiIiIBEbhUkREREQCo3ApIiIiIoFRuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQmMwqWIiIiIBEbhUkREREQCo3ApIiIiIoFRuBQRERGRwChcioiIiEhgZgyXZnadmfWY2SNF59rM7E4z684dW3Pnzcy+ZmZPmNlDZvbKhWy8iIiIiNSWanourwdOKzv3KeDfnHNdwL/lXgOcDnTlHh8GvhFMM0VERERkOZgxXDrn7gb6yk6fAdyQe34D8KdF5290WfcBLWa2PqjGioiIiEhtm+ucy3XOuRcAcsfO3PkNwHNF1+3InRMRERGRVSAc8P2swjk33Qe6u7sDbkJtfJeIiIjIYlqsnNPV1TXt+3MNl7vMbL1z7oXcsHdP7vwOYFPRdRuB5+fTwKB0d3cv2neJiIiILKZayjlzHRbfBpybe34u8LOi8+fkqsaPBgbzw+ciIiIisvLN2HNpZjcDJwFrzWwH8FngSuCHZnY+8CzwrtzltwNvAp4ARoHzFqDNIiIiIlKjZgyXzrmzpnjrlArXOuDC+TZKRERERJYn7dAjIiIiIoFRuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQmMwqWIiIiIBEbhUkREREQCo3ApIiIiIoFRuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISmPBSN0BEREREZikxivX14vX14hqbl7o1JRQuRURERGpJUXC0vl6sryf3vAfrz50fHSlcnjrhTXDC25ewwaUULkVEREQWy/gYtqcHr78X29ObDYt7sqGxECKLguNUXCSCa+vEb+3A32fTIjS8egqXIiIiIkGYNjj24vX1VB8cWzvw2zpxbR3ZEJk7urYO/LYOaGgGs4kPdXcv4A+bHYVLERERkZmMj00ExL7eomHrnmCDY2sHNJYFx2VG4VJERERWt5mCY38vtnd4xtvkg2O2d7EzGyLbs0fX3rkigmM1FC5FRERk5SoEx+JexuyQdWEIW8ExUAqXIiIisjwtRHBsLZ7bOHFUcKyewqWIiIjUnvFEaUFMPiz2FRXHVBMcw7ng2K7guFjmFS7N7BPABwEHPAycB6wHfgC0AQ8C73POJefZThEREVkpZgyOvdjeoRlvMzk4Ti6QcY0tCo6LbM7h0sw2AB8HDnXOjZnZD4EzgTcBVzvnfmBm3wTOB74RSGtFRESktuWDY3kl9VyDY27pHQXHUs5PY15tDkDPt1VhoM7MUkA98AJwMnB27v0bgM+hcCkiIrL8LUZwzBXIuIZm8LxF+FG1w/kZXGoAl+zPPsb7cck+vHgH4X1OAcAfe5GxBy7Aoi3UH3PdEre4sjmHS+fcTjP7MvAsMAb8X+C3wIBzLp27bAewYd6tFBERkYVVEhwrbDm4p8rgGAqX9TJ2TKqsXk3B0TkH/jgWimdfZxKkd94+ESCTfbhkP/54P6QGyc40LOW1HlkIlxZphMwoLhVazJ8xK+bc5B9R1QfNWoFbgf8BDAA/yr3+rHPugNw1m4DbnXMvL/7s4OBg4Uu7a2hFeRERkZXIUuNEhweIDPYRHe4nMpR9RIf6iAz3Ex3qJzy2d8b7+F6IVFMryaY2Uo0tE88L51pJr2kAW/nB0fxxPH8YLzNEKDOUPfrDpCIbSNRvASA6/kfae64hGdvMns6Lcp9Lsn7nX1a8p8PwvQb8UCOZUBO+10Qm1EQ6sp6xNUflLnKYP4rz6pdsSkBXV00Nn/IAABfSSURBVFfheXNz86RGzGdY/FTgKedcL4CZ3QYcC7SYWTjXe7kReL7aBi6k7u7uRfsuERGRRZMcn1joe09P5b2qR2bR41jcy1gYti6a4+h5GBDNPdYs9O9bRM5Pg5/EwvXZ18l+Us/fkRui7ivqbeyHzFjFe4T3eT2xrncB4O+NMtaTpi5amneSoXdAuBGLtmGxViyae0RaMG9uPZK1lHPmEy6fBY42s3qyw+KnAL8Bfgm8k2zF+LnAz+bbSBERkVVpUYJjLjzmguNK45wPqeGSIej8MHSo7ZWE2/8EgPSuuxj/w5WE1r2O+GGXZD+bHiP15A2Vb+xFcqGwbSIcRlvxmg4sXGL1G6k/8aeFIfG86AEfWpgfWyPmM+fyfjP7MdnlhtLA74BvA/8b+IGZfSF37togGioiIrKiFAfH8jmO+UKZaoNj69rJe1Sv8ODo0mOAj4Wzfaf+6E7SL/68UART0svoMhXvYV4EcuHSIk2AB/7E6okWayOy77tLexdzYZLwGmyGYWmzEIRqd27kQpnznMv5KJ5zuVhqqbtYRERWuOT41MUxQQbH1g5cU+uKCY7OT+GSA6XhcDx7DO9zCqHmgwFIPvMjUn+8lsi+7yR6wAcByPQ/ROJ3f135xuGGkt5Fi2UDYqj5MEIth+W+OwOWC4TL0FLlnKDnXIqIiKw+UwbH3PP+Xmx4cMbbTATHid1iXFt2j2rXvnKCY3ZYeig7jJzrZcwM/TfpXXeVVkyP90N66h13vDX7FcKlRZvBi+KKeiStfgOR/d83KUBapAULRWds51znOspkCpciIiJ5k4Jj8TqOCo55zrnscjhFPYvFvY2Rfd+Ft2YTAMnub5LesY1o10eIbHp79vOjO0k/d9vkG5uXDYPF4TD3CLUcXrgsvM+phPd5fcmwtBdrJ7r/exb2h0tVFC5FRGR1qBAciwtkvL6e2QfH/LqNrcUhcvkGR+cncckBLNxQqJjO9D1IuvfXhd7FwjxGf3zK+4TWHlMIl9n5iY3g/ML7XlMXkZedP1EEE8vOZSTSWNWwtK2CpY6WM4VLERFZ/pLjWP/ukt1i5hYcQxNzHFdIcHQuA6kh/AqFLi7ZT6zrz7FoCwDjf7iSTO+viR3+acKdrwUgM/xH0jv/ZfKNvVhR72IbFm0pPPca9i9cFtnvTKKbzyr9aP1Govu9a+F+tCwphUsREaltUwXH4gKZ4YEZbzMpOFbaq3oZBUeX3osb78PiHYWlbtIv/pJM34Nlu78MAv6U9/E3vYNQLlxarAOLrS2prg61HYmFLphULW3huqraOVNFtaw8CpciIrJ0UslJvYyFParnEhxby+Y4LrPg6DLjk3oWixfvjh3+mULhSeL3f4M/9DjxP/kKoeZDAcgMPU76xTsn3zjSVFotHW3LDUW34sU7CpfFDvwoHPjRko+GGg8g1HjAwv1oWXEULkVEZGFMGxxzBTLzDY75oeumFqjhal9/vA+X7MNbs38hHKae+ymZwUdLq6Uzo9Pex6UGsVgbAFa3HksNg58uvB9edwJew+aJHsZYKxZpzq7nKLJIFC5FRGT2Ukmsf3eF9RvnEBxb1lK+zaDf1lnTwdE5B+mRir2LLtmP81PED7+0cH3iNx/Hje+m7tgbsXgnAJmBR8j03lt6YwsXLaUzefeXfJENQPywT01qV6j50EIvpshSUbgUEZFSMwXH/l68of4Zb+M8r+L+1IXg2NaBa26tqeDonMMlXsQlBwtrKgIkn7wBf+SpkqV3cKmpb2Qezl1SqGr21uyLizRCZqLCOrLhLYQ7jisJkoQbNEdRlj2FSxGR1aQQHHMV1Lk9q7N7V88hOLZ24LdPDE9nC2U6ayo4Oj+NSw0UbQs4ULIuo4UbiB18Ue5qn7H/OB/wqT/pXzAv+7/JTN/v8Ye2l944VF+5dzE3HA0Tm9HFt3xxUrtCbVsW5geLLDGFSxGRlWK64JhfEHwFBUfnp3GjO3CZREkv4/j2q/ETPRPrMaamX4LIYhMFLWYhvMaXgXmQGQOvEYDI/u8Bf7x02Z1chbaIlFK4FBFZDlJJbGBPUVgs2qN6rsExv25j2S4ySxkcXSZRNn+xr2Qo2luzqbCXtEsNMvbAn2PRVuqPv7lwj0zf73DjPUV3taI1GEt7F7PP20vaUPfqr09qV7j9VQvye0VWIoVLEZGlNlNw7O/FG6wyOLasLVqCp3aCo8sk8Pc+CzhCTQcVzo0/+uWSoWoyY9Pfp6gX0iItWP0mLNaOc64wVzF60AWlhTGRZu0bLbKIFC5FRBZSPjiW71GdL5SZd3DMhcbWDlxL26IFR+d8SA0Vehj98mrpZD+h1i1EN58JgL/3GRK/uQiv8QDqXn1N9iZelMzu/yhZsBsvOjGHMVY0BJ3vbcxVWgOYF6L+6O9Malt47dEL+ttFZHoKlyIic5VOlc5xnHdwLN1mcMmCY2oYf+RpCMUINR0IgJ/oJfn41yaGqFP9JXtFV2LhhonnsXa8hpfhrdlv4px52UXBww2FIWpC9aqWFlnmFC5FRCopCY5FYbGwGHjP/INjvkBmgYOj85MTFdKTehizi3eH159KZMObAcgMPMz4w58n1P4aQq/4ewDMC5PZ85+lNw43lPUyllZOe3X7FC71YmupO+p/TWpbuOOYBfvdIrI0FC5FZPWZMTj2YkP9mHPT3saZh2ttn9irur2ztMI6P8cxtHD/qfXH9+CPPI1FWwk1vhSAzPAfSXZ/uxAeSY/MfJ9cDyWAxTvxmg/FW7PvxAWRJmJHfK4oQDZjXjTw3yMiy5/CpYisLPng2N+Ll1u3sXwx8GqDo9+2dprg2IFrbgs0ODrnIDM6Ze9i/nVk85mEO18LQKb330n+9z8R3vBmQgd9LHcnH3/gvyZubKGySun8DjBF6zPWrS9cHmo8gLo/+UpJ28xCmssoIlVRuBSR5SPo4Fhpr+q2hQmOAP7o8/h7n8ar31joFUzv+S2pp74/sSajn5z5PmMvFp579RvxWrfg1W0oOreJ+JYvFoaoiTQWdooREVloCpciUhvS6WxY7C/en7q0UKaWgqPzM9ldXwp7S/dP7m1M9hM76CJCrS/P/sQX/g+pZ24hsv85RPc/O3ejFP7QYxM3DtXlehNbSnd+KdoJxivuZWx7JXVtryxpm4XihMrOiYgsFoVLEVl46TQ2ULxX9fyDY3GBTHGhTFA9jpnhP+L2PovXfEihMCX1ws9JP3cr/nh+15fp2wuULObtNbyUUPtrSsNh82HEX/n/TYTIcN282y4ispQULkVkfgrBsWy3mNzWg9bfiw32VRccc8UxQQdHl0nkqqX7ynoZ+0p6HuNbLsdbsyn7s577CekXf0704E9OVD374/gjT+XuahBpLhS4eLGy3V/yvY3xdYV2hNedSHjdiSVts0gjoZbDZv2bRERq1bzCpZm1AN8FDif7V/gPAI8DtwCbgaeBdzvnZl6vQ0Rqz0IFx5ICmaJ1HOfY45gZeAR/dCeh9lfjxdoASD7zI9Iv3IEb74fMaFX3cck+yIVLr/kQQv44FpvYGjDccSxe04HZ8Bhpxjz9/VxEpNx8/8v4VeAO59w7zSwK1AN/A/ybc+5KM/sU8Cngknl+j4gErTg4FoXFiWHrnlkGx44KldWzC47OOUgPV6yW9sf7S3of64+5DguvASD51E34/b8n9orLC+GSzChudGf2uUXK1mIs2/klXzUdW1toS2TDmwvrPuZZtJVQtHUWf8giIqvPnMOlmTUBJwDvB3DOJYGkmZ0BnJS77AbgLhQuRRZXOo0NFu1VHWRwLCmQ6Zx1j2N69/24xC7C+7y+ML9w/L+/Qab333HJAXDpqu7jkgOFcBlqeyVebC0WaSq8H97wZsLrXpetlg43aNcXEZFFYm6G/7lM+UGzLcC3gUeBVwC/BS4CdjrnWoqu63fOlfxVf3BwsPCl3d3dc/p+kVUrkyYyMkh0qJ/IUD/RoT4iw/nn/dnnw4PYDMUmDiPV2EyqqZVkYyupprbs86ZWUo2tJJvaSDU0TR0cXRovM0zIH8bLDBHKDOP5Q7nnuaM/jPnj7NpweeFjHS98gUh6Fz37/A3pSLawpaVvK/V77wPAtzoyoUb8UBO+10Qm1IQfasweC6+b8L0G0PI6IiKLrqurq/C8ubl50t/c5xMuXwXcBxznnLvfzL4KDAEfm024XCzd3d0lfxgiNSmTxgb2TLFXde4x2IfNsKezM8M1t+NyC35XXJKnuR3CpcHROVfo4XPOkdn1C1xygPCmPyucTzx8GZmBhyE1VPXPqj/xp1goDkDyyRtwqSEi+76rUCjjJ3rB+dmh6ZB2fRERma2lyjmVwuV85lzuAHY45+7Pvf4x2fmVu8xsvXPuBTNbD/RMeQeR1aQkOE4s+u3lth60PdUHRz+3V3Vhj+ryyuqi4Jjd9WWsdP3F0T/iBiZXSxNppP413wTAzBh//H9BZpTw+jdApDH3O8ZywdLLbgE4xVqMxbvA4MUK7Y++9NxJv8mLdwTzZywiIktuzuHSOfeimT1nZgc55x4HTiE7RP4ocC5wZe74s0BaKlLLMmlsoG/SbjELGRydy2AWyn4ukyD94i/Af4JI++GF+43958fw9z4L/nh1v6PsuvBL3jjpkujBn8hWSUeaCt8vIiKSN99q8Y8BW3OV4k8C5wEe8EMzOx94FnjXPL9DZGmVB8dcgYzXX7Tl4MAcg2OusjpfYe03tYIbLaqOzj+6cckHcM9M9D56jQdQd+SVuZv7JB//GngxwhvPmBja9pPZwOjFsj2KsQo7v5TsMd1S0uZY10cm/Q71MoqIyHTmFS6dc78HXlXhrVPmc1+RRZMPjv29ePl1G+ccHIsXAJ8IjpmWtbi2dmhdl+1xTA6Q3nUXeJHCUjfOTzP263NwqQGY4bsK35kcmHgRqiP8ktOxSHP287kexfiWy7FQfXZLQVVLi4jIItAKwLJyTRscc+FxLsExv/VgayuZ5jr8hhB+1OEyg2X7TP8x+/yFfsKZE4h1XJy9X2qYZPc3sboNhXBpXjjby+j87HBztHQ9xkq7vxTmQJKdHxk7+KJJbfeK1m0UERFZDAqXsjxl0tk5jGV7VJfsVT2b4JjvZWxdi2trg7b1+G2dZNY4UuN/wOrWEln/eiBb2Tz2wEdhbATGqmuuSw0XnlusnfCGt+LVrSu5pu4138YijZgXmd2fhYiISA1RuJTaM1Nw7O/F+vdUHRz91rX4He1k2hvJtMTxGyL4dYYfzeC8BC41WDSv8TdE9nsn0ZdlZ3b4fb8j9fgNeC1HFMKlRRohPQIWqryXdHRyL2N+sXAAC9cTO+jCSe0t7CwjIiKyjClcyuLyM4XleKZckmdgD+bPHBwzra24tg5oWYff1kmq3SO1Zg/W9FLCG0/FtbST2ftHEr/7K/CfmfhwMveY6t6ZROG51W8ksu87sDX7TZwLxal/7Q9zu75oEW8REZFiCpcSnOmCY75AZobg6AwydUZmbQv+2iYyrfVkmmL49R4u5uOHkviM4TJDkBkleuBbiGx8GwCpF35OcvuXCUVihNZmF+e2SCP4KQjVl63FWLYeY359xkhzdpmdHC/eQfSAD01qZ/E2gyIiIjJB4VKq42cKy/FkC2QqLMkzRXB0gIuBpQFn+M1tjO/fSPIlRtjbRHjNYbjWDsbrXiTRd1PuE4nco0gm98izMC4zsS5jqPEAIpvPxmuc2KHA6vYp2R1GREREFpbCpVQRHHuxgd2TgqML53oZ4+A3GH6HkWltINMUx18Two+DH0nhLAH4xPf9GKHNb4BwhOTTPyD15PVE9n0t7oAzszfs/y/oAyItRdXRldZjzK7TmB2Wnlhex2vYTLRhc0kbzUIQ0kLfIiIii0XhcqXLB8f+yXtUFwpkioKjM/DrwEuA5bJkYl+P5AEesd5GwuyDa+1gdN+9jDVvr/CFaWBk8unwGlxDPYSzldCh1i3w0vcTapnYTcZrPpz6k/4F8xQGRUREliuFy+WsyuCI7+Ni4NdZtqexLvvcbzH8fQy/LkRmTQy/znDRbKJsSLyZUOuh+G0dJMbuINX3CzJvO4/Ipj/NfvWLv4DHnijqWcz3KBbv9lL0CMVKmh5qPphQ88El5xQqRURElj+Fy1rlZ7DB/tJK6qJCGfp2wcgeXMwRGnGYy35stCtEaq2xZneG0ED25NBxDYwdkJ7hC13u4WHRFtJHnohrPQKAUJ+PtR+I13JE4erQuhOpX/c67foiIiIiJRQul8I0wdH6e3CjPbhUPy7qF3oaU3WW7W1sN/yNkIkbRKIAtP7HS/AaXpLtZdznEVKh5+Do8whvOBnXupb0zp/CM7dMsZd08S4wbRBpzM5TLBJq20KobUvJufJrREREREDhMniF4Jgbqu7vhT27YOgFGO4h8vxgdo5jJsPIy0NkWjwafpMilNvpZeCECOP7h6jqH40XxaJtJP7i7/HWbMqe2n0/0fE9WNuRuLr1AET2fSfR/d61QD9YREREZILC5WyUBUfr34kb2IHbuwvGduOnB3HsJRNz2V7G/KMDWGfgOzofGccc+E2tjB+QId2UJNJ4HKGGA3Ftnfje/dj4w1i8vah3sbxiOnskVD9pWDq89jWTmq2haxEREVksCpd5ueDInueh7xls4AXCfZlCocze9c/ihxI0/ypZqKLuOy1Kap0Ha4tvVHm42Ihj0Wb2XvG3sHZfiEQJ9dyLl0mQOf5V+NEWAKKcSnRhf6mIiIjIglk14TIytBO3/RkYfA6GX8CN7cZP9uMywzgbw4+k8ePgYtlePm/c0fGvEwt0j2/J7hKT7mjGi6/DtXdgDS/g+cNYqDHby1jfCY3rsbq1Zb2NLZgXmdSmcOfxi/b7RURERBbDqgmX657/EqNrciXVIaCh/Irc0LEPXiaChetJvOut0LYOv62DiPcUNDSTOOFoLFwPQCT3EBEREZGsVRMuSdYTckk8F8e8xmxvYl0HNOwDLZuw1v2wNZ25amkPyC4HnhdmS+X7ioiIiEjBqgmXO464kq6urpkvFBEREZE585a6ASIiIiKycihcioiIiEhgFC5FREREJDAKlyIiIiISmHmHSzMLmdnvzOxfcq/3N7P7zazbzG4xM60JLiIiIrJKBNFzeRGwvej1l4CrnXNdQD9wfgDfISIiIiLLwLzCpZltBN4MfDf32oCTgR/nLrkB+NP5fIeIiIiILB/zXefyH4G/Bhpzr9uBAedcfv3xHcCG6W7Q3d09zyZUbzG/S0RERGQxLVbOmWnd8DmHSzN7C9DjnPutmZ2UP13hUjfdfRZrYfPu7m4toi4iIiIrUi3lnPn0XB4HvM3M3gTEgSayPZktZhbO9V5uBJ6ffzPnr1b+wEVERESCVks5Z85zLp1zlzrnNjrnNgNnAr9wzr0H+CXwztxl5wI/m3crRURERGRZWIh1Li8BPmlmT5Cdg3ntAnyHiIiIiNQgc27aKZEiIiIiIlXTDj0iIiIiEpgVES7NzJnZ94teh82sN79rkIiIiMhKYWYjM7x/l5m9arHaU25FhEtgL3C4mdXlXr8e2LmE7RERERFZlVZKuAT4V7K7BQGcBdycf8PMjjKzX+f2QP+1mR2UO3+PmW0puu7fzeyIRW21iIiIyCyZ2UnFI7Rmdo2ZvX8Jm1SwksLlD4AzzSwOHAHcX/TeY8AJzrkjgb8Dvpg7/13g/QBmdiAQc849tGgtFhEREVlhVky4zIXCzWR7LW8ve7sZ+JGZPQJcDRyWO/8j4C1mFgE+AFy/KI0VERERWaFWTLjM2QZ8maIh8ZzLgF865w4H3kp2RyGcc6PAncAZwLuB/3/xmioiIiIyZ2lKc1x8qRpSbj7bP9ai64BB59zDRfudQ7bnMl/g8/6yz3wX+GfgHudc34K3UERERGT+ngEONbMY2WB5CnDv0jYpa0X1XDrndjjnvlrhrauAK8zs34FQ2Wd+CwwB31uEJoqIiIjMmZmFgXHn3HPAD4GHgK3A75a0YUVW/Q49ZvYS4C7gYOecv8TNEREREZmSmb0C+I5z7qilbstUVlTP5WyZ2Tlkq8o/rWApIiIitczM/pxsXclnlrot01n1PZciIiIiEpxV3XMpIiIiIsFSuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQnM/wPzU/wQ6wjIYgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "plt.plot([\"May\", \"Jul\"], [fl_before, fl_after], label=\"FL\", lw=2)\n",
    "plt.plot([\"May\", \"Jul\"], [poa_before, poa_after], label=\"POA\", lw=2)\n",
    "\n",
    "plt.plot([\"May\", \"Jul\"], [poa_before, poa_before+(fl_after-fl_before)],\n",
    "         label=\"Counterfactual\", lw=2, color=\"C2\", ls=\"-.\")\n",
    "\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See that small difference between the red and the yellow dashed lines? If you really focus you can see the small treatment effect on Porto Alegre. \n",
    "\n",
    "![img](./data/img/diff-in-diff/cant-read.png)\n",
    "\n",
    "\n",
    "Now, what you might be asking yourself is \"how much can I trust this estimator? It is my right to have standard errors reported to me!\". Which makes sense, since estimators without them look silly. To do so, we will use a neat trick that uses regression. Specifically, we will estimate the following linear model\n",
    "\n",
    "$\n",
    "Y_i = \\beta_0 + \\beta_1 POA_i + \\beta_2 Jul_i + \\beta_3 POA_i*Jul_i + e_i\n",
    "$\n",
    "\n",
    "Notice that \\\\(\\beta_0\\\\) is the baseline of the control. In our case, is the level of deposits in Florianopolis in the month of May. If we turn on the treated city dummy, we get \\\\(\\beta_1\\\\). So \\\\(\\beta_0 + \\beta_1\\\\) is the baseline of Porto Alegre in May, before the intervention, and \\\\(\\beta_1\\\\) is the increase of Porto Alegre baseline on top of Florianopolis. If we turn the POA dummy off and turn the July Dummy on, we get \\\\(\\beta_0 + \\beta_2\\\\), which is the level of Florianópolis in July, after the intervention period. \\\\(\\beta_2\\\\) is then the trend of the control, since we add it on top of the baseline to get the level of the control at the period post intervention. As a recap, \\\\(\\beta_1\\\\) is the increment we get by going from the treated to the control, \\\\(\\beta_2\\\\) is the increment we get by going from the period before to the period after the intervention. Finally, if we turn both dummies on, we get \\\\(\\beta_3\\\\). \\\\(\\beta_0 + \\beta_1 + \\beta_2 + \\beta_3\\\\) is the level in Porto Alegre after the intervention. So \\\\(\\beta_3\\\\) is the incremental impact when you go from May to July and from Florianopolis to POA. In other words, it is the Difference in Difference estimator. \n",
    "\n",
    "If you don't believe me, check for yourself. You should get the exact same number we got above. And also notice how we get our much wanted standard errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>  171.6423</td> <td>    2.363</td> <td>   72.625</td> <td> 0.000</td> <td>  167.009</td> <td>  176.276</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>poa</th>       <td> -125.6263</td> <td>    4.484</td> <td>  -28.015</td> <td> 0.000</td> <td> -134.418</td> <td> -116.835</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>jul</th>       <td>   34.5232</td> <td>    3.036</td> <td>   11.372</td> <td> 0.000</td> <td>   28.571</td> <td>   40.475</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>poa:jul</th>   <td>    6.5246</td> <td>    5.729</td> <td>    1.139</td> <td> 0.255</td> <td>   -4.706</td> <td>   17.755</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smf.ols('deposits ~ poa*jul', data=data).fit().summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non Parallel Trends\n",
    "\n",
    "One obvious problem with Diff-in-Diff is failure to satisfy the parallel trend assumption. If the growth trend from the treated is different from the trend of the control, diff-in-diff will be biased. This is a common problem with non-random data, where the decision to treat a region is based on its potential to respond well to the treatment, or when the treatment is targeted at regions that are not performing very well. Take our marketing example. We decided to test billboards in Porto Alegre not in order to check the effect of billboards in general. The reason is simply because sales perform poorly there. Perhaps online marketing is not working there. In this case, it could be that the growth we would see in Porto Alegre without a billboard would be lower than the growth we observe in other cities. This would cause us to underestimate the effect of the billboard there.\n",
    "\n",
    "One way to check if this is happening is to plot the trend using past periods. For example, let's suppose POA had a small decreasing trend but Florianopolis was on a steep ascent. In this case, showing periods from before would reveal those trends and we would know Diff-in-Diff is not a reliable estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApcAAAErCAYAAACLjKgCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxc9X3v/9d3du2bF9mSsVlk8IJZvIAXwNjOQpqUmybtDdnI1rRJmuampZdwk5T8SlMS2ps0Db8maQIFGhpICL2Q24QGG4zBNgZsCIsNiMUYyZZk2dql2b/3jzMzmpFGliyNRiPp/Xw89NDonDMz3zHH1pvv9jHWWkREREREcsE11Q0QERERkZlD4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREcsYzFW/a1dWlJeoiIiIi01xFRYUZekw9lyIiIiKSMwqXIiIiIpIzsyZcNjY2TnUTRPJO973MRrrvZTYqpPt+1oRLEREREZl8CpciIiIikjNTslp8JNZaent7icfjOX/tQCBAV1dXzl93tnO5XJSWlmLMsMViIiIiMgsVVLjs7e3F7/fj8/ly/tp+v59AIJDz153twuEwvb29lJWVTXVTREREpAAU1LB4PB6flGApk8fn801KT7OIiIhMTwUVLkVERERkeiuoYfFCUF1dzfLly1M/33333Rw5coRbb72Ve++9dwpbJiIiIlL4FC6HKCoq4oknnsg4duTIkSlqjYiIiMjIrLVECmx2msKliIiIyDTSHY6z61iIR5pDbG8O8oEzi/hwxVS3alBBh8vKf23O6et1frJu1GsGBgbYtGkTAIsXL+buu+/OaRtEREREToe1lhdORtiRCJP7WsNE7eD5p4+HFS4LWbZhcREREZF8OhGM8ejRENubgjx6NETrwODYt8vAJfN8bK3zs60uwAU1Xt54vXMKW5upoMPlWHoaxyoYDObstURERERyKRq3HGgPs705xI6mIAfaI6R1TrKw2MXWugDb6gNcscBPpb9wN/wp6HApIiIiMlMd7YuxoznIjuYQO48G6QwPxkmfCzbU+tm60M/W+gDLKj3TphqewuUYPfbYYxlbFN1xxx2sW7duClskIiIi00koZnmyNZSaO3mwI5px/qwyN1vrA2yrC7Cp1keJt3B7J09F4XKI5ubhi4guu+wyWlpapqA1IiIiMp292R1le3OQ7c0hHj8Woj9tJU6Jx3DZAj9b6/xsrQtwVvnMiGWjfgpjzCLgLqAWiAP/Yq39njGmGrgXWAIcBv7IWtthnD7b7wHvAfqBT1hrD0xO80VEREQKR18kzuMtTu/kjqYgb/TEMs4vr/KwrS7A1roAl8734XdPj6Hu0zGWiBwF/tJae8AYUwbsN8Y8DHwC2GGt/ZYx5ivAV4DrgauAhsTXJcAPEt9FREREZhRrLYc6o+xocnon97aGCKdtal7hM2xZGGBrvZ8tCwMsLHFPXWPzZNRwaa09BhxLPO4xxhwC6oCrgc2Jy+4EduKEy6uBu6y1FnjSGFNpjFmQeB0RERGRaa0zFGfnUWfe5CPNQY72D6ZJA6ye403MnfRz8RwfHtfM6508ldMa3DfGLAEuAvYB85OB0Vp7zBgzL3FZHfB22tOaEseyhsvGxsbU40AggN/vP50mnRZtRzQ5uru7aWtrm+pmyAjS/46JzBa67yWX4hYO9brY2+Fib4ebF3tcxBkMjNVey/qqGBuqYqyrjFHpBeiCLnizK3/tzNd939DQcMrzYw6XxphS4JfA/7DWdp9iOXy2EzbLsWEN7OrqIhAIjLVJpyUYDE7aa8925eXlLFq0aKqbIVk0NjaO+o+AyEyj+15yobU/xiNHQ+xoDvJoc4gTocHeSY+BDfN9bKsLsKXOz8pqL64p3iaokO77MYVLY4wXJ1jeba29P3G4NTncbYxZACS7rpqA9KRRDxzNVYNFREREci0St+xrC/NIc5DtTSGePxnJOH9GqTuxEMfPZQv8lPum5zZB+TCW1eIGuA04ZK39TtqpB4FrgW8lvj+QdvzPjDH34Czk6ZpO8y2rq6tZvnw5sViMpUuX8oMf/IDi4mKam5u57rrreOWVV4jH47zrXe/ipptuwufzpZ57zTXX0N7ezsMPPzyFn0BERETG4q2eKI8k9pzcdSxET2RwoDXghstqnQ3Mt9b5Oad8+mxiPtXGErs3Ah8Dthhjnkt8vQcnVL7DGNMIvCPxM8CvgTeA14AfA5/PfbMnT7K2+N69e/H5fNx+++1Ya/nYxz7G7/3e73HgwAH2799PX18fN910U+p5nZ2dPP/883R1dXH48OGp+wAiIiKS1UDUsr0pyA37Oll3fysX3NfKl/d28p9HgvRELOdVevj8ihLuf2cNb354Ib945xz+dHkpDRVeBcvTMJbV4k+QfR4lwNYs11vgCxNsV0FYv349L730Ert27cLv9/PRj34UALfbzd/93d9xwQUXcMMNN1BcXMyDDz7Iu9/9bubOncv999/PX/zFX0xx60VERGY3ay2NXVGnXndzkN0tIYJp206Wew1XLPSn5k4uKp0Zm5hPtYL+Uyy9dnPuXgvovXPnmK+PRqM8/PDDbNu2jUOHDnHhhRdmnC8vL6e+vp433niDlStX8stf/pLrr7+eefPm8fGPf1zhUkREZAp0h+M8dszZwHzH0RBv92ZuYn5hjTdVEWftPB/eWbZNUD4UdLicCgMDA2zatAlwei4/9rGPcdttt2XtDrfWYoyhra2NN954g/Xr12OMwePxcPDgwYxa5CIiIpJ7cWt54WTEqdfdFOSptjBpFRap8bucMFkf4MqFfuYVzfxNzKdaQYfL0+lpHE0wGGQsGxEl51ymW7ZsGQ8++GDGse7ubpqbmznzzDO566676OrqYtWqVQD09PRw//33K1yKiIhMgvZgjEebk5uYhzgeHNwmyG1g/XwfWxMruy+omfptgmabgg6XheKKK67gG9/4Bj/72c+45ppriMVifO1rX+PDH/4wxcXF/PKXv+S+++5j3bp1ABw+fJj3v//9fO1rX5vilouIiEx/0bhl//Fwau7ks+2RjA2064rdbK13hrqvWOCn0q9tgqaSwuUYGGP46U9/yl/+5V/y93//98Tjcd7xjnfw13/917z11ls0NTWxdu3a1PVLliyhrKyMZ555hjVr1kxhy0VERKan5r4YO5qD7GgOsvNoiK7wYJz0uWBjrT81d/K8Sm0TVEgULodobm7Oery+vp5777132PHFixdz6NChYcd37dqV87aJiIjMVKGYZW9riO1NIR5pDnKwM5px/pxyD1vqnJXdG2t9lHjVO1moFC5FRERkSrzRHWV7k9M7+XhLmP60lTglHsPlC/xsSwx3LylTZJku9F9KRERE8qI3EufxY6FUVZw3ezK3CVpZ7WXrQmdl96XzfPjcGuqejhQuRUREZFJYaznYEWVHc5DtzSH2toaIDC7sptJn2JJY1b2lLsCCYm0TNBMoXIqIiEjOdITi7DzqhMlHmoMc6x9MkwZYO9fL1roA2+oDXFTjxa1NzGcchUsREREZt1jc8tyJCNubg+xoCvFMe5h42j5B84tcTpis87N5oZ/qgHonZzqFSxERETktrf3JbYJCPHI0SEdoME16XbBxvo9t9QG21gVYUaVtgmYbhcshWltbueGGGzhw4AB+v58zzjiDm2++mXPOOScnr//444/j8/m45JJLTvu5n/70pzl06BAf+chH+MIXvjDm53V2dnLffffxmc985rTfM+lzn/sc7373u7n66qvH/RoiIjI9hWOWfW1hHknMnXzhZCTj/OJSN++oD7Clzs9lC/yUaZugWU3hMo21lo9+9KNcc8013H777QA8//zztLW15SxcPvHEE5SWlp5WuIxGo5w4cYJ9+/bx4osvnvZ7dnV1cdttt00oXIqIyOxyuCeaWtW962iI3rRtgorchssW+BLD3QHOKnerd1JSFC7T7Nq1C4/Hw6c+9anUsVWrVmGt5etf/zrbt2/HGMN1113HH/zBH/D4449z6623pjZX/6u/+isuvPBCPvKRj3D++edzzTXX8NBDDxGNRrnjjjvw+/3867/+K263m3vvvZdbbrmFpUuX8uUvf5mmpiYAbr75Zi699FJuvvlmWlpaOHLkCDU1NRw8eJD29nY2bdrELbfcQmNjI3fccQfhcJizzjqLH/3oRxQXF9PW1saXv/xlDh8+DMB3vvMdfvSjH/Hmm2+yadMmrrzySt75zneO2O5vf/vbPPTQQwSDQdatW8c//uM/6h8MEZFZoD8aZ3dLmO1NQR45GqKxK3MT82WVnlS97vXz/QQ8+t0g2RV0uOx75N2ndb2r7ByK1t467PklWx4a0/MPHTrEhRdeOOz4gw8+yAsvvMATTzzBiRMn2LJlCxs2bBj19Wpqati1axc/+clP+P73v8/3v/99PvnJT1JaWsoXv/hFAD7zmc/w+c9/nvXr1/P222/zgQ98gKeeegqA5557joceeoiioiLeeustPvShD/HEE08AcN5553HttdcC8Ld/+7f827/9G3/yJ3/C9ddfz8aNG7n77ruJxWL09vZy4403cujQodRzH3/88RHb/NnPfpbrr78+9fihhx7iqquuGtOfn4iITB/WWl7tijr1upuC7G4NEUrbdrLcZ9i8wM+2+gBbFvqpLy3oyCAFRHfKGDz55JN84AMfwO12M2/ePDZs2MCBAwcoKys75fPe9773AXDhhRfyq1/9Kus1O3fu5OWXX0793NPTQ09PDwBXXXUVRUVFWZ938OBBvvnNb9LV1UVvby9bt24FnN7XH/7whwC43W4qKiro7Owc82fdtWsX//RP/8TAwAAdHR0sW7ZM4VJEZIboCsd57GgotRinqS9zE/OL5nhTK7vXzPXh0TZBMg4FHS7H2uOYq+cvW7aMBx54YNhxa22Wq8Hj8RCPD+7fFQwGM877/X7ACXnRaObwQlI8Hufhhx/OGiJLSkpGbOvnP/957r77bs4//3zuvvvuVK/kWIzU7mAwyHXXXcejjz5KfX09N99887DPJCIi00fcWp4/EWFHYu7kU21hYmm/0uYEXKl63Vcu9DO3SNsEycRpOVeayy+/nHA4zJ133pk6duDAASorK/mP//gPYrEY7e3t7Nmzh9WrV7No0SJefvllQqEQXV1dPPbYY6O+R2lpaapnEmDLli38+Mc/Tv38/PPPj6mtvb291NbWEolE+MUvfpE6fsUVV3DbbbcBEIvF6O7upqysLOM9R2p3MkjW1NTQ29vLgw8+OKa2iIhI4WgPxvj56/18dtdJlt7TwuZfHeemA93sbQ0DsH6+j69fXM7O983l1Q/V8i+XV/NHZxcrWErOFHTPZb4ZY/jpT3/KDTfcwHe/+10CgUBqK6K+vj42bdqEMYa/+Zu/Yf78+QC8//3vZ+PGjZx99tmsWrVq1Pe46qqr+PjHP86vf/1rbrnlFr797W9z3XXXsWHDBmKxGBs2bOC73/3uqK/z1a9+la1bt7Jo0SKWL19Ob28vAN/61rf40pe+xE9/+lNcLhff+c53WLduHZdeeinr169n27Zt3HTTTVnbXVlZybXXXsuGDRs444wzuOiiiybwpykiIvkQjVuePh5mR7Mz3P1ce4T08bb6EjfbEuUVr1jop8KnfiWZXGakId/J1NXVlfVNu7q6qKiomJT3DAaDBAKBSXnt2W4y/7vJxDQ2NtLQ0DDVzRDJq9lw3zf1RnnkaIjtTUF2HgvRHR78tep3w8b5frbWO3Mnl1ZoE/PZYKru+4qKimE3l3ouRUREClwwatnbGkrV6z7UmTmPv6HCw9Y6P1vrAmys9VHsUe+kTJ1Rw6Ux5nbgvUCbtXZl4tiFwA+BABAFPm+tfco4/2v0PeA9QD/wCWvtgclqvIiIyExkreWN7phTr7s5yOPHwgykrcQp9RguX+gsxNlS52dJmfqKpHCM5W68A7gVuCvt2C3A/2et/Y0x5j2JnzcDVwENia9LgB8kvouIiMgp9ETiPH4slJo7ebgnc5ug86u9bKtzhrvXzfXhc2uoWwrTqOHSWrvLGLNk6GGgPPG4AjiaeHw1cJd1JnI+aYypNMYssNYey1F7RUREZgRrLS92RJ163U1BnmwLExncJY5qv7NN0JaFzmKc2mKt5pbpYbz96P8D+C9jzD/gbGeULFdTB7yddl1T4tiI4bKxsXGwMR5n0rHP5xtns05NezbmXjgc5uTJk7S1tU11U2QE6X/HRGaLQr3vuyKwr9PN3g43T3a6aA8Pzo10YTm/LM76qhgbquKcVxrHbZydQHqaoWekFxVJyNd9P9rCofGGy88BX7bW/tIY80fAbcA2IFsf/SmXo6c30FpLb28voVBonM0aWXd3N+Xl5aNfKKfF5XKxZMkSrUQsULNh1azIUIV038filgPtkURFnCD72yPE034rLih2JSriONsEVfm1EEfGp5Du+/GGy2uBLyUe/wL4SeJxE7Ao7bp6BofMR2WMGbWk4ni1tbWxaNGi0S8UERGZgJb+WKq84qNHg3SEBtOk1wWbav3O3Mm6AMurtE2QzDzjDZdHgSuAncAWINkP+yDwZ8aYe3AW8nRpvqWIiMxk4ZjlybYwO5qC7Dga4sWTkYzzS8rcvKMuwNZ6P5tq/ZR61TspM9tYtiL6Gc5K8DnGmCbgRuCPge8ZYzxAEPhs4vJf42xD9BrOVkSfnIQ2i4iITKnDPVF2NAfZ3hTi8WMheqODvZPFHsNltT5nuLs+wFnl2iZIZpexrBa/ZoRTq7Nca4EvTLRRIiIihaQ/GueJY+HUvpOvd2duE7S80pOqiHPpfD9+bRMks5j+d0pERGQIay2vdEXZ3uTMndzTGiKUlicrfIYrFzobmG+tC1BXom2CRJIULkVERIDOUJzHjjkbmO9oCtHcP5gmDXDxHG9iZbef1XN9eFzqnRTJRuFSRERmpbi1/O5EhO1NQR45GuKptjBpFRaZG3Cxtc7PtvoAVy70UxNQ76TIWChciojIrHF8IMYjR0PsSATK9uBgSRyPgQ3zfWyrD7C1zs/51V5c2iZI5LQpXIqIyIwViVuebgvzSHOI7c1BnjuRuU3QolJ3as/Jyxf4KfdpmyCRiVK4FBGRGaUlaNjzSh/bm4M8djREd2RwrDvgdjYx35KYO9lQoU3MRXJN4VJERKaluLUc7onxUkeEg4mvF09GeL27COhMXbe0wpOaO7lhvp8ij8KkyGRSuBQRkYJ3fCDGwY4IL3VEU0Hy5c4o/WmblyeVuC1X1hWxtc7ZKmhxmX7VieST/saJiEjB6I/GebkjmtYb6YTJ42kLb9ItKHaxvMqb9uXB1/4Wy86tz3PLRSRJ4VJERPIuFre80RPlYDJInnTC5Js9MYb3RUKZ17Cs0gmPy6u8LK/2sqLKS5V/+AKcxpOT334RGZnCpYiITBprLW0D8cSQ9uCw9iudEYKx4dd7DDRUeDJ6IpdXeTmj1K2FNyLThMKliIjkRG8kzsudTnh86eTgsPaJUPYh7foS92BPZOKrocKjutwi05zCpYiInJZo3PJ6d3TYApvDPVm6IoFyn2HFkJ7IZZVeKrMMaYvI9KdwKSIiWVlrOdYfT4XHlxI9ka92RQhlyZFelzOkvXLIApu6Eg1pi8wmCpciIkJ3OM6htNXZydXaneFsy2ucyjbLq7ysSBvWbqjw4HUpRIrMdgqXIiKzSCRuaeyKpvVGOo/f7s0+pF3pM4kQOdgTuazKqzKJIjIihUsRkRnIWktzXyzVE5nsjXy1K0oky/oanwvOTWz1syJtWHtBsUtD2iJyWhQuRUSmuc5QnEOdmZuOv9QRoXuEIe0lZe5UeEwOa59d7sGjIW0RyQGFSxGRaSIcs7yaNqSdDJNNfdmHtGv8rtTq7BXVTpg8r9JDqVdD2iIyeRQuRUQKjLWWI72xjJ7Igx0RGruiZCmlTcAN51UOzolMDmvPK9KQtojkn8KliMgU6gjFM8ofHuyIcqgzQk9keIo0wFnJIe3q5CIbD2eVeXBrSFtECoTCpYhIHoRillc6I8MW2Bzrz169Zm7AlbHp+IoqL+dWeijRkLaIFLhRw6Ux5nbgvUCbtXZl2vEvAn8GRIH/tNb+z8TxG4BPAzHgz621/zUZDRcRKUTxxJB2evnDgx0RXuuOEssypF3sMZxX6Rm2wGZukTv/jRcRyYGx9FzeAdwK3JU8YIy5ErgaWGWtDRlj5iWOLwc+BKwAFgLbjTFLrbXZZ5uLiExjJ4KxjPKHBzsiHOqI0pdlYqTLONVr0mtpr6jysqTMjUvzIkVkBhk1XFprdxljlgw5/DngW9baUOKatsTxq4F7EsffNMa8BqwD9uasxSIieTYQdYa0XxqywKZ1IPuQdm2RK6P84fIqL+dWeinyKESKyMw33jmXS4HLjDHfBILAddbap4E64Mm065oSx0bU2Ng4ziacvny+l0ih0H0/djELzUHDa30uXu93vr/W76JpwBBneDAsclnOLolzTrHlnJI45xTHObskTqV3yIUd0NSRn88gDt33Mhvl675vaGg45fnxhksPUAVcCqwFfm6MOQuy/OsL2XfxHWMDc6WxsTFv7yVSKHTfj+z4gLPVz4tpPZEvd0QZyDIx0m3g3HJP5gKbai9nlGpIuxDpvpfZqJDu+/GGyybgfmutBZ4yxsSBOYnji9KuqweOTqyJIiLj1x+N83JHNDGkPTisfTyYfUh7YXH6kLYTJpdWeAloSFtEZEzGGy7/D7AF2GmMWQr4gHbgQeDfjTHfwVnQ0wA8lYuGioicSixueaMnysFkkEys1n6zJ5Z1+KTMazJ6IpNfVX5t9SMiMhFj2YroZ8BmYI4xpgm4EbgduN0Y8yIQBq5N9GK+ZIz5OXAQZ4uiL2iluIjkkrWW1oF4ap/IZE/kK50Rgln+tfEYWFrhYXl15gKbRSVuVa8REZkEY1ktfs0Ipz46wvXfBL45kUaJiAD0RuIcSoTH9GHtk6HsQ9r1Je7UPpHJr4YKDz63QqSISL6oQo+ITLlo3PJadzQxlD04P/Kt3uwDH+U+k6qfnSyBuKzKS4VPQ9oiIlNN4VJE8sZay7H+eEb5w5c6orzaGSGcpTPS63KGtFcMWWBTpyFtEZGCpXApIpOiOxznUNqcyGRvZGc4++5kZ5S6M8ofLq/yck6FB69LIVJEZDpRuBSRCYnELY1d0bTeSOfx2yMMaVf5TUb5w+VVHs6r9FKuIW0RkRlB4VJExsRaS1NfLKP84UsdERq7okSyDGn73XBuhRMeV1R5U6u1a4tcGtIWEZnBFC5FZJjO0OC8yFSY7IzQnWVI2wBnlrmH9UaeVe7BoyFtEZFZR+FSZBYLxSyvJoe0Tw6Gyeb+7EPaNX5XRvnD5VVezqv0UOrVkLaIiDgULkVmAWstR3pjmT2RiSHtaJb1NQE3nFc5uDo7uVp7noa0RURkFAqXIjOMtU5v5J6WMLve9NL8ynEOdUboiWQf0j673J2x6fiKKg9nlnlwa0hbRETGQeFSZJqLW8tLHVH2tITY3RJiT2uY9mByhY0Xp0IrzCtyZZQ/XFHl5dxKD8UeDWmLiEjuKFyKTDPRuOX5ExF2t4TY3Rpmb2uIriELbeYXudhY62cxXWxeupAV1R7mBNxT1GIREZlNFC5FClwoZnm2PczuljB7WkPsaw3TO2Si5KJSNxvm+9hY62fjfD9nlTsVbBobT9Cw0D9FLRcRkdlI4VKkwPRH4zzdFmFPqzPM/czxMMEhi7fPLnezYb6fjbV+NtT6OKNUf5VFRKQw6DeSyBTrDsd5qi2cmi95oD08bFPyZZUeJ0jO97Gh1k9tsYa4RUSkMClciuRZRyjO3tZQapj7dycixNNGuV0GLqjxpoa518/3UaP5kiIiMk0oXIpMsraBGHtawuxODHMf7IhmnPcYWD3Xm+iZ9HPJfB8VqrMtIiLTlMKlSI419UbZ0zo4zN3YlRkm/W5YPSex+KbWx9q5PkpU4UZERGYIhUuRCbDWcrgnluiVDLOnJcRbvZmrb4o9hkvm+VLzJVfP8RHwaINyERGZmRQuRU5DsvpNcr7knpYQR/szV9+Uew3rE0FyY62fC2q8eFXtRkREZgmFS5FTiFvLiycj7Gl1eiUzq984qv2u1OKbDbU+VlZ5VTpRRERmLYVLkTTRuOV3JyLsaQnxRGuYJ09R/WZjrY8N8/2cW+nBZRQmRUREQOFSZrlQzHIgWf2mJcS+tjB9WarfbEwOc6dVvxEREZHhRg2XxpjbgfcCbdbalUPOXQf8PTDXWttunN+43wPeA/QDn7DWHsh9s0XGx6l+E2Z3Ypg7W/Wbc8o9bEj0Sqr6jYiIyOkZy2/NO4BbgbvSDxpjFgHvAI6kHb4KaEh8XQL8IPFdZEp0h+Psa3OC5O6WMM+eGF79ZnmlJ7H4xsf6+ap+IyIiMhGjhktr7S5jzJIsp74L/E/ggbRjVwN3WWst8KQxptIYs8BaeywXjRUZzclgjL2t4dQ+k8+fzF79JjlfUtVvREREcmtc433GmN8Hmq21vxsy96wOeDvt56bEsRHDZWNj43iaMC75fC/Jj/YwPNft5kCXi2e73LzWn7kZudtYzi+Lc1FFnIvLY1xQHic1yh2Gk2/Dyfw3O69038tspPteZqN83fcNDQ2nPH/a4dIYUwx8FXhnttNZjtksx1JGa2CuNDY25u29ZPI09UZT8yVHqn6zZq7TK7mp1seaWV79Rve9zEa672U2KqT7fjw9l2cDZwLJXst64IAxZh1OT+WitGvrgaMTbaTMTtZa3uyJpcoo7m4JcWRI9ZsSj2FdovrNxlo/F6v6jYiIyJQ67XBprX0BmJf82RhzGFiTWC3+IPBnxph7cBbydGm+pYyVtZZXuqLsaUnW5Q5xbGj1G59h/bzkhuWqfiMiIlJoxrIV0c+AzcAcY0wTcKO19rYRLv81zjZEr+FsRfTJHLVTZqBY3PJSRyStlGKYE6HMMFnjd2VsC6TqNyIiIoVtLKvFrxnl/JK0xxb4wsSbJTNRJG55/kSE3S0hdreG2dsaontI9ZvaRPWbDbVO7+TSClW/ERERmU60O7RMmlDMsv/44LZAT2WpfnNGqTs1X3JjrZ8zy1T9RkREZMA20voAAB78SURBVDpTuJSc6YvEeea4U/1md6L6TWiE6jcba/1smO9jkarfiIiIzCj6zS7jlqx+s7vFmS95oD3MkI5Jlld6UsPcG+b7ma/qNyIiIjOawqWM2clgjD2tzuKb3S1hXshS/ebCGu/gApz5PqpV/UZERGRWUbiUEbX2x1KruHe3hDjYmblhucckNixPDHOvm+ejwjd7NywXERERhUtJ83ZvNLX4Zk9LmNe6s1e/2VjrZ+N8Vb8RERGR4RQuZ6lk9ZsnWkLsSWwN9HaW6jeXzPOxIbH4ZvVcH363VnKLiIjIyBQuZ4lk9Ztkr+SI1W/mO72SG2v9rFL1GxERETlNCpczVCxuebEjkpovubf11NVvNtb6WKHqNyIiIjJBCpczRCRu+d2JiDPE3RJib1t4xOo3ya2Bzq3waMNyERERySmFy2kqGLUcaE8svmkNj1j9JrlZuarfiIiISD4oXE4TfZE4Tx8PszsxzL2/fXj1m4YKDxvmDy7AUfUbERERyTeljwLVFY6zL7ktUGuIZ9sjw6vfVHnYON8Z5l4/36fqNyIiIjLlFC4LRLL6TXKY+1TVbzbOd8Kkqt+IiIhIoVG4nCIt/TH2JILknlNUv9lY6wxzXzLPR7mq34iIiEiBU7jMk7d7o+xuSdblDvF6d+aEyYA7WUrR2Wdy7TwfxR6FSREREZleFC4ngbWWN7pj7E4EyT2jVL/ZWOvj4jmqfiMiIiLTn8JlDsSt5ZXOaCpI7mkJ0TKQvfrNpsRqblW/ERERkZlI4XIcktVvdreEU/MmT45Q/Sa5z6Sq34iIiMhsoHA5BpG45bn2SGq+5JOtYbojmfsCLSh2JYKkM8y9VNVvREREZBZSuMwiGLXsb3d6JXcnqt/0D9lkcnGpOzVfcuN8P0tU/UZERERk9HBpjLkdeC/QZq1dmTj298D7gDDwOvBJa21n4twNwKeBGPDn1tr/mqS250yy+s0TiWHuU1W/SQ5z16v6jYiIiMgwY0lIdwC3AnelHXsYuMFaGzXGfBu4AbjeGLMc+BCwAlgIbDfGLLXWDolqU6srHOfJ1uR8yVNUv6n1pzYsV/UbERERkdGNGi6ttbuMMUuGHPtt2o9PAh9MPL4auMdaGwLeNMa8BqwD9uakteMUi1sebXdz+8lOdreEebEje/WbZK/khlo/VX7tMSkiIiJyunIxtvsp4N7E4zqcsJnUlDg2pVwGvv26jxORPgC8Llg715ca5l6n6jciIiIiOTGhcGmM+SoQBe5OHspymc1yLKWxsXEiTRizP1zgIWoNF1fEOL8sTsDd75wYgNa3oDUvrRDJv3z9HRMpJLrvZTbK133f0NBwyvPjDpfGmGtxFvpstdYmA2QTsCjtsnrg6EQamCufpjFv7yVSKBobdd/L7KP7XmajQrrvxzUWbIx5N3A98PvW2v60Uw8CHzLG+I0xZwINwFMTb6aIiIiITAdj2YroZ8BmYI4xpgm4EWd1uB94OLG345PW2j+11r5kjPk5cBBnuPwLhbZSXEREREQmz1hWi1+T5fBtp7j+m8A3J9IoEREREZmetERaRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJm1HBpjLndGNNmjHkx7Vi1MeZhY0xj4ntV4rgxxvyTMeY1Y8zzxpiLJ7PxIiIiIlJYxtJzeQfw7iHHvgLssNY2ADsSPwNcBTQkvj4L/CA3zRQRERGR6WDUcGmt3QWcHHL4auDOxOM7gf+Wdvwu63gSqDTGLMhVY0VERESksI13zuV8a+0xgMT3eYnjdcDbadc1JY6JiIiIyCzgyfHrmSzH7Kme0NjYmOMmFMZ7iRQK3fcyG+m+l9koX/d9Q0PDKc+PN1y2GmMWWGuPJYa92xLHm4BFadfVA0cn0sBcaWxszNt7iRQK3fcyG+m+l9mokO778Q6LPwhcm3h8LfBA2vGPJ1aNXwp0JYfPRURERGTmG7Xn0hjzM2AzMMcY0wTcCHwL+Lkx5tPAEeAPE5f/GngP8BrQD3xyEtosIiIiIgVq1HBprb1mhFNbs1xrgS9MtFEiIiIiMj2pQo+IiIiI5IzCpYiIiIjkjMKliIiIiOSMwqWIiIiI5IzCpYiIiIjkjMKliIiIiOSMwqWIiIiI5IzCpYiIiIjkjMKliIiIiOSMwqWIiIiI5IzCpYiIiIjkjMKliIiIiOSMwqWIiIiI5IxnqhsgIiIiIuPU3YmJhqe6FRkULkVERESmk+5OPPt34XlqJ+5DzxHd/F7Y+L6pblWKwqWIiIhIoRsSKI2NA2DdbgiHprhxmRQuRURERApRdyee/Y/jeepR3C8/h4kPBsroikuIrt1M9OKNUFoOjY1T3NhBCpciIiIihSIZKJ/eifvQs5mBctWQQFmgFC5FREREptIMCJTpFC5FRERE8u1UgfL8dUTXXTmtAmU6hUsRERGRfOjpxLP/CWcOZdZAuZnoxZumZaBMp3ApIiIiMllSgXIn7kMHZmygTDehcGmM+TLwGcACLwCfBBYA9wDVwAHgY9bawtrdU0RERGSyzMJAmW7c4dIYUwf8ObDcWjtgjPk58CHgPcB3rbX3GGN+CHwa+EFOWisiIiJSiE4ZKNcm5lDmPlDaSC82Xlh9eBMdFvcARcaYCFAMHAO2AB9OnL8T+AYKlyIiIjLT9HbheebxvAZKa+PEe14jduIZYif3E+8+hKfuvWDekbP3mKhxh0trbbMx5h+AI8AA8FtgP9BprY0mLmsC6ibcShEREZFCkAyUTz+G++D+4YFy7WaiqzdBaUXO3zp8+F4ib98Pka7Bg8YN0T7w5vztxm0iw+JVwNXAmUAn8AvgqiyX2lO9TmMed5TP53uJFArd9zIb6b6XXHL391L5yrNUHnqGsjdfHiy9aFx0n7WCjmWr6Tr3ImLFpc4TjrUBbRN6T2/4CIH+3zFQfBFRXz0Apd0nKY90EXVXEQosJxRYRiiwFOsqAvJ33zc0NJzy/ESGxbcBb1prjwMYY+4HNgCVxhhPoveyHjg6kQbmSmNjY97eS6RQ6L6X2Uj3veTESD2ULhfRFWudRTmrN+EqraAGqJng28WDbRhPCcZTAkDold8S7fkt1TVz8Z11pfPe4TnYyHsxxfUYYzKeX0j3/UTC5RHgUmNMMc6w+FbgGeBR4IM4K8avBR6YaCNFREREJl1v1+CinKGBcuVgoMzFkLeNhYl1vkDs5DPETuzH9h/Bd95f4F34TgA88y4Dlwd3zZrUc4yvCuOrmvB7T7aJzLncZ4y5D2e7oSjwLPAvwH8C9xhj/jZx7LZcNFREREQk55KB8umduA8ewMRiQFqgXHuFEyjLKif0NtZabH8TsZP7ncU4nS9APDR4gbsYoj2DP1atwl21akLvOVUmtFrcWnsjcOOQw28A6ybyuiIiIiKTJl+BMtpPrON3id7JZ7DB1ozzrtKzcdeswV29BlfFMoxrZtS2mRmfQkRERORU8hQok2y0j/7HPwQ2MnjQW467ejXu6tV4alZPiyHu8VC4FBERkZmptxvPgUQt76GBcsWawTmUE+6h7CP86g+I979NYPU/YozBeEpwlZ0NxuUEypo1uMrOwRh3Lj5ZQVO4FBERkZkjFSgTi3JyHChtPEa85xXivW/hrUvswOguInriGYh0YvuPYEoWAxC4+H9jXDM/TA6lcCkiIiLT2yQHyniondiJ/c7cyZPPQrQXjAvPvMsw3lKMceFf9mVcgfmY4jNSz5uNwRIULkVERGQ6msRAaeNh4p0vETv5DNET+7F9hzPOm6KFuGvWYOMhDM7G6Z45l0z4I80UCpciIiIyPfT1JPahfDRLoFzt1PIeZ6CMB48Ta9/r9FB2PDdkm6AA7qoLcFevwV29Glfxwlx9ohlJ4VJEREQKVzJQPr0T90vPDA+UazcTXX0ZlJ9eoLSxIMRCGJ+zIXqsfR/hV/85dd5VemZiZfcaXJXLMS5f7j7TDKdwKSIiIoUlI1Dux8SiwMQDZVKk+TeEX/1nPPW/j7/hjwGc/SbnXe70TtZcjMs/J2cfZ7ZRuBQREZGpN1KgNOMPlDbSQ6zjWWInnsFVuQrvgm0AzrC2jWLDJ1PXuopqCaz8X7n9TLOUwqWIiIhMjb6ewUU5QwPl8osTcyjHHiitjRPvaXTKK57cT7zrZcCpD+4OdwyGy4oVFG+6JzUkLrmlcCkiIiL5kwyUTz+G+8VnsgTKzcRWX4YtH1v1mnjoJLGTBwa3CYp0DZ40blwVq3DXrMFTs2bwsMsDCpaTRuFSREREJleOAyVA+I1/I9b+JPHe1zOOm8D8VL1ud9UFGE9xTj+KjE7hUkRERHIvh4EyPtBCrONZPAvemSqfGOt8wQmWLl9imyCnZrcprscYM6kfTU5N4VJERERyo68Hz7O7nTmUEwiUNh7O2Pon+OxXsMEWXKVn4S4/FwDfkg9j7X/HXbkS4/ZP3meS06ZwKSIiIuOXg0BprcX2NznzJk88Q6zzRYrW347LXwOAZ97lxAeaAVfqOe7qCyf1Y8n4KVyKiIjI6RktUK69gtiay08dKKN9xDqeS9XstsG2jPPx7ldwzd0AgO+cT03eZ5GcU7gUERGR0U0wUFobJ977Rto2QQfBxgYv8Fbgrr4YT80a3NUXY3xjX9wjhUXhUkRERLLr68Hz7B6nlvfQQLnsosEh74rqU75MrPMlQi/+LTbcMXjQuHBVrHAW4tSswVV2Dsa4Rn4RmTYULkVERGRQMlA+vRP3C0+fdqCMD7QSPfZf4C7Ct/gPAXAV12HDHRj/nER5xTW4qy7EeEvz9rEkfxQuRUREZrvRAuXazcTWZA+U8eBxbKQHd9lZznMiXUQO/zsmMA/vGR/EGIPxVVJ06e2YogXaJmgWULgUERGZjfp78RzY7QTKF5/BRCPA6IHSxsPEO18kemI/sZP7sX2HcVWspGj1PwDgKjsH7xkfwFV5AWABJ0y6ihfm89PJFFK4FBERmS36ezPnUI4hUFprsQNHUwtxYh2/g3ho8DXdAYyvAmvjGOPCGBe+c/44359MCsiEwqUxphL4CbAS539PPgW8AtwLLAEOA39kre0Y4SVERERkMqUC5U7cLz6dGSjPu5DouiuHB8roALHO3zmB8sR+bPBYxku6Ss9KzJ1cjatiOcblzetHksI20Z7L7wEPWWs/aIzxAcXA/wJ2WGu/ZYz5CvAV4PoJvo+IiIiM1YiB0iQCZWJRTqWzSbm1FhI9jwChg7cQa987+Hqe0sSqbqfEYnJzc5Fsxh0ujTHlwOXAJwCstWEgbIy5GticuOxOYCcFEC7n7fkNvmcfxZZXYssqsaUV2PIK53FZJfgDU91EERGR8TvNQJkUfvNuos2/xr/8OtzVFwHgrlmDDXc4q7qr1+Aqb0jV9BYZjbHWju+JxlwI/AtwELgA2A98CWi21lamXddhrc3YCbWrqyv1po2NjeN6/9O17AdfJ3CiZcTzMa+PaHEZseJSosVlRFPf0x+XEi1xjsX8RaAVbyIiMoVcoQEqXv0dVQefoeyNl3AlV3lj6F28lM5lq+k872KipRVg43jDR/AHDzJQcgkxjxMyyzofoKxnOz3l76an4vem8uPINNHQ0JB6XFFRMSwMTSRcrgGeBDZaa/cZY74HdANfPJ1wmS/tv7yLBV6D6enE9HRhujsxvV2JnzsxkchpvZ51u53ez7LKRG9oxZCfK6GswjleVoktLQe31k9JfjU2Nmb8IyAyG8z4+36gL7HK+zHcLzyV0UMZP3cVkXVXpnoo46ETxE4eSCzGOQDRHgB8Sz+Pt/73AYgPtECsH1NyprYJmsam6r7PFi4nknaagCZr7b7Ez/fhzK9sNcYssNYeM8YsANpGfIU86li1njkj/aFbC8GBwaDZ05X5PT2Idie+B/sxXSeh6+SY3t8aA8VlzlB86eBwfCqYlqV/dx7j8+fwT0BERKatUwTK2HkXEF27meiay4mXlxPvOkis/QFir+4n3vt6xsuYQG2iGs7S1DFXUW1eP4rMfOMOl9baFmPM28aYc621rwBbcYbIDwLXAt9KfH8gJy2dTMZAUTG2qBg7b4z7cIVDmN7uIWF0SDDt7kwdp68H09eN6esG3h7TW9hA0fDe0SwhNBVGi0o0VC8iMlOkB8oXn0qNsA0LlKVFRFsfJXbkVmIdz0FsYPA1XH7cVatSVXFM0UL1Tsqkm+g47ReBuxMrxd8APgm4gJ8bYz4NHAH+cILvUZh8fmz1XGz13LFdH49Bb09aCE0G0bTe0Z5ED2mqd3QAExyA9pHniqazHm/mQqXSisEFTFmCKaVl4NIEbRGRgjHQl7YoJ3ugjFx8CTFvMFURh2gf4Vf/GWwMAFOyGHf1Gjw1q3FVrMS4fVP1aWSWmlC4tNY+B6zJcmrrRF53RnK5oTzRAzmW662Fgb7hITRjmD6zd9SEgpjOduhsH1OTrDFQWj4kiGbvFU31jnr1j5SISE4lA+XTO50h7/RAee4FRNc5PZS2soZ48DgDT34K3EUUb/oZxrgxnhK8i/87JjDXWdkdGGOnh8gk0QqTQmUMFJdii0ux8+vH9pxwaGwhNNk72tcNieA6VjZQnBimP1UITYTo0goIaFW9iMgwowTKyLr1hM6uJBp6BRt+jkDl+wEw/jkY3xyMtxQb7sQk9pv0nfXxKfsoIkMpXM4kPj+2Zh62Zt7Yro9FM+aNMnSIvqcrEUTTQmmwHxPsh+NHx/QW1utNC6HJAFqR/Vh5JRSXgcs1gT8EEZECdYpAGT13FaG1Kwgv8hLtO0i8+y54M556ajx0Epe/GmMMRZf8EOPWgk8pXAqXs5nbg62ozij5dUrWQn/v8N7RoSE0fWFTOIQ5eRxOHh/bW7hczhB9aQWUjxBCM3pHy8GjsmMiUqAG+vE8l6jlPSRQRlasIHjxGYTnhIn2vAiRVyH5/+3GhatiZWIT89UYX2qHPwVLKXgKlzJ2xkBJGbakDFs7xqH6UDAtiHalDcuPsO1Tfy+muwO6Owb/kR2FLS4ZvXc0uQVUeSX4i8b/ZyAiMppUoNyJ+4V9mUPeS1fRv+4MgmWvE+9/A3gdEjvaGf/cVJh0V1+E8ZRM3WcQmQCFS5lc/gDWX4udM8Z91KLRwfDZ25XZO9rTOTh0n77Qqb8P098Hrc1jegvr848+RJ++GX5xqeaNisipjRAoo8WG0CWLYcl6WPNBbNUcose2Ez+0HVxe3JXnD24TVLxI2wTJjKBwKYXF48FW1gyrfTuieDwxVD+W3tHOwaH6E61wonVMb+FUYyofFkIZOkSvakwis0uWQGldEPcYWLqK6LrN9NefIHzsl3jqIvir5gDgmbMOc8FNuCvPx7gDU/whRHJPvwFlenO5nO2USsuxC84Y/XprITSQuYI+FURH2Ax/oA/T1QFdHWNuli0pG3klfWLxUvo+pKrGJDJNDPTjeW5vYg7lPohEiJUb+s9xE2qoIVIZxLfwA3iXfQoAV08j7shR3JXnp17CeMvx1Kydqk8gMukULmV2MQYCxc6WSnMXjO05kfDgoqX03tFsc0Z7OqG3G9PXg+nrgZYxVmPyB04dQof0jqoak0geJQPl0ztxP78Pa8OEF7joX+0mtLiMeCCSuLAPgDiD27u5yxpwr7pxChotMnUULkVG4/WdfjWmvp7h+4sO3Qw//VgoiAm1jL0ak9szvFc0I4g634uOn8CU+JyeUa8P6008dnsUTkVOJS1Qup7fR6w0QqjORWiLm8i8gFOLDoAIeMpwV1+cWIxzMS7/GKf1iMxQCpciueZyQ7L3ceHi0a+3FoL92UuAjhRMgwOYzhPQeeKUL33eSG9pjFNtyevDetODp2/wcSKIJh/bREAd9pzE87JfNyTYev3gdivYSmEa0kNpImEscOL9PmLl6VNXXLjKz3XCZM0aXGXnYIxK6YokKVyKTDVjoKgEW1SCnV83tueEQ9mrLw0JoaGuDgLGQCTsDO+HQ873WBTCIed16JnczzeENa7B8OkbKYQOBthhgfVUAdg79Dp/RgDWQisZJpicQ+n0UPYvj9E/30VVLEJs6flE126G6qcxwbdx16x2AmXVRRhv2VS3XKRg6V9akenI58dWz8NWn7oaU2NjIw0NDcNPxGMQiUAkhAkngmckNBhCI2EIpz8OOd+HhtRIGCLpj8Op10u/Lv31TSwG4SCEg5i+SfrzGYF1uYYF1lSPqneE3tdEmB3+nMzwOrzHNjMA41LPVsFIBErXsw8Taz9A4LXw4KmGcmIlYbq/8TXci68AwB95F3hKtE2QyBgpXIrMRi43+N3OPqSJQ/aUT8ihWBSikYzAmhFIh4XZtMfhtJA6WgDOuC7kPDceh1DQ2dw/X583wbrdI0xFGBpsfWO8boQAnAy9accUbIFgP+5nn4CDvyY6cIi+Wkv0XBec68IdOBd7wTuJrrkcT+xVPMaNq+rC1FONt3QKGy4y/Shcikh+uT3Ol79oaoJtKqRm610dGmZPHYCzPmdoz244cT4Wg9gABAemINh6sofUIVMMsk1FGGsATr92sGfX52wXNlWC/bgO/Ib4G78lEj9M/3ywZxog8WU9uCtWEvzcn+AqPRMAD2NcuCciI1K4FJHZIxlsA8X5DbbWQiw2fPpBePjjzF7aLFMWwiMH4PQQnDFlIRaFWBQT7M/Hp8386B7vCGE2M4Rm9sQOnX4wUo/t0IVlfqyNUP3KfxB77iUi3hZi5QYWAolI77JVuOddirvuclyVKzAuX97/TERmOoVLEZHJZgx4PE4FqqLi1OH8BdvoqcNsRiDNDK/DenZHCLPZpiKYSBgTjUA0ghnIzwTb4BkuAlf6CBYDGEzUhce1GNfiLbgWbcYVUM+kyGRTuBQRmcmMAY8XPF5sUUnqcN6C7ZD5tRmLv4ZNL8gegIc9Z+iUhbTX8XaGcHW58JatxHXuezF1l2I051QkrxQuRURkcqTvp5p2eLKDbdNIuySISF5M4UxrEREREZlpFC5FREREJGcULkVEREQkZxQuRURERCRnJhwujTFuY8yzxpj/m/j5TGPMPmNMozHmXmOMNhETERERmSVy0XP5JeBQ2s/fBr5rrW0AOoBP5+A9RERERGQamFC4NMbUA78H/CTxswG2APclLrkT+G8TeQ8RERERmT6MtePfccwYcx9wM1AGXAd8AnjSWntO4vwi4DfW2pXpz+vq6kq9aWNj47jfX0RERETyK30f2YqKCjP0/Lg3UTfGvBdos9buN8ZsTh7Ocukp02u+Nrpt1Ka6MgvpvpfZSPe9zEaFdN+Pu+fSGHMz8DEgCgSAcuA/gHcBtdbaqDFmPfANa+270p+b3nMpIiIiItNTtp7Lcc+5tNbeYK2tt9YuAT4EPGKt/QjwKPDBxGXXAg+M9z1EREREZHqZjH0urwf+whjzGlAD3DYJ7yEiIiIiBWhCC3pERERERNKpQo+IiIiI5MyMC5fGmN6pboNIvhhjrDHm39J+9hhjjicrZonMJLrfRRyjZR1jzE5jzJp8tWeoGRcuRWaZPmClMaYo8fM7gObTeQFjzLi3JBPJswnf7yIy+WZkuDTGlBpjdhhjDhhjXjDGXJ04vsQYc8gY82NjzEvGmN+m/SMlMl39BqdSFsA1wM+SJ4wx64wxe4wxzya+n5s4/gljzC+MMb8Cfpv/JouM23ju98eNMRemXbfbGLMqr60WyTFjzOb0XntjzK3GmE9MYZNSZmS4BILA+621FwNXAv87UZoSoAH4/621K4BO4ANT1EaRXLkH+JAxJgCsAvalnXsZuNxaexHw18DfpZ1bD1xrrd2St5aKTNx47vef4FSQwxizFPBba5/PW4tFZpmZOhxmgL8zxlwOxIE6YH7i3JvW2ucSj/cDS/LfPJHcsdY+b4xZgtOL8+shpyuAO40xDTjVsrxp5x621p7MSyNFcmSc9/svgK8bY/4K+BRwR14aKzJLzdSey48Ac4HV1toLgVacKkIAobTrYszcgC2zy4PAP5A2RJhwE/CotXYl8D4G/x6AM39NZDo6rfvdWtsPPAxcDfwR8O/5a6rIpImSmeMCI12YbzM1WFXg1D2PGGOuBBZPdYNEJtntQJe19gVjzOa04xUMLnj4RL4bJTJJxnO//wT4FfC4euxlhngLWG6M8eMEy63AE1PbJMeM6rlMrHoNAXcDa4wxz+D0Yr48pQ0TmWTW2iZr7feynLoFuNkYsxtw57lZIpNiPPe7tXY/0A38ax6aKDJpklnHWvs28HPgeZzc8+yUNizNjKrQY4y5APixtXbdVLdFREQKhzFmIbATOM9aG5/i5oiM23TIOjOm59IY86c482++NtVtERGRwmGM+TjOqvKvKljKdDZdss6M6rkUERERkak1Y3ouRURERGTqKVyKiIiISM4oXIqIiIhIzihcioiIiEjOKFyKiIiISM4oXIqIiIhIzvw/TYTpMouBtxgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "x = [\"Jan\", \"Mar\", \"May\", \"Jul\"]\n",
    "\n",
    "plt.plot(x, [120, 150, fl_before,  fl_after], label=\"FL\", lw=2)\n",
    "plt.plot(x, [60, 50, poa_before, poa_after], label=\"POA\", lw=2)\n",
    "\n",
    "plt.plot([\"May\", \"Jul\"], [poa_before, poa_before+(fl_after-fl_before)], label=\"Counterfactual\", lw=2, color=\"C2\", ls=\"-.\")\n",
    "\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will see how to solve this problem with synthetic control. It will use multiple cities to create a synthetic city that closely follows the trend of the city of interest. But for now, remember that you always need to check if you have parallel trends when applying diff-in-diff. \n",
    "\n",
    "![img](./data/img/diff-in-diff/non-parallel.png)\n",
    "\n",
    "One final issue that it's worth mentioning is that you won't be able to place confidence intervals around your Diff-in-Diff estimator if you only have aggregated thata. Say for instance you don't have data on what each of our customers from Florianópolis or Porto Alegre did. Instead, you only have the average deposits before and after the intervention for both cities. In this case, you will still be able to estimate the causal effect by Diff-in-Diff, but you won't know the variance of it. That's because all the variability in your data got squashed out in aggregation.\n",
    "\n",
    "## Key Ideas\n",
    "\n",
    "We've explored a technique widely applied when we are estimating causal effects at more macro entities (schools, cities, states, countries...). Difference in Difference takes a treated unit before and after the treatment and compares the trend in the outcome to that of a control unit. Here, we've seen how this could be applied at estimating the effect of a city specific marketing camping.\n",
    "\n",
    "Finally, we looked at how Diff-in-Diff fails if the trend between the treated and control unit is not the same. We also saw how diff-in-diff will be problematic if we only have aggregated data.\n",
    "\n",
    "## References\n",
    "\n",
    "I like to think of this entire book as a tribute to Joshua Angrist, Alberto Abadie and Christopher Walters for their amazing Econometrics class. Most of the ideas here are taken from their classes at the American Economic Association. Watching them is what is keeping me sane during this tough year of 2020.\n",
    "* [Cross-Section Econometrics](https://www.aeaweb.org/conference/cont-ed/2017-webcasts)\n",
    "* [Mastering Mostly Harmless Econometrics](https://www.aeaweb.org/conference/cont-ed/2020-webcasts)\n",
    "\n",
    "I'll also like to reference the amazing books from Angrist. They have shown me that Econometrics, or 'Metrics as they call it, is not only extremely useful but also profoundly fun.\n",
    "\n",
    "* [Mostly Harmless Econometrics](https://www.mostlyharmlesseconometrics.com/)\n",
    "* [Mastering 'Metrics](https://www.masteringmetrics.com/)\n",
    "\n",
    "Other important reference is Miguel Hernan and Jamie Robins' book. It has been my trustworthy companion in the most thorny causal questions I had to answer.\n",
    "\n",
    "* [Causal Inference Book](https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/)\n",
    "\n",
    "Finally, I'd also like to compliment Scott Cunningham and his brilliant work mingling Causal Inference and Rap quotes:\n",
    "\n",
    "* [Causal Inference: The Mixtape](https://www.scunning.com/mixtape.html)\n",
    "\n",
    "![img](./data/img/poetry.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
